Joan Solà

Quaternion kinematics for the error-state Kalman filter

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1 Quaternion definition and properties

Definition of quaternion

One introduction to the quaternion that I find particularly attractive is given by the Cayley-Dickson construction: If we have two complex numbers A = a + bi and C = c + di, then constructing Q = A + Cj and defining k ij yields a number in the space of quaternions H,

Q = a + bi + cj + dk ∈ H , (1)
where {a, b, c, d} ∈ R, and {i, j, k} are three imaginary unit numbers defined so that

i 2 = j 2 = k 2 = ijk = -1 , (2a)
from which we can derive

ij = -ji = k , jk = -kj = i , ki = -ik = j . (2b)
From ( 1) we see that we can embed complex numbers, and thus real and imaginary numbers, in the quaternion definition, in the sense that real, imaginary and complex numbers are indeed quaternions,

Q = a ∈ R ⊂ H , Q = bi ∈ I ⊂ H , Q = a + bi ∈ Z ⊂ H . (3)
Likewise, and for the sake of completeness, we may define numbers in the tri-dimensional imaginary subspace of H. We refer to them as pure quaternions, and may note H p = Im(H) the space of pure quaternions,

Q = bi + cj + dk ∈ H p ⊂ H . (4)
It is noticeable that, while regular complex numbers of unit length z = e iθ can encode rotations in the 2D plane (with one complex product, x = z • x), "extended complex numbers" or quaternions of unit length q = e (uxi+uyj+uzk)θ/2 encode rotations in the 3D space (with a double quaternion product, x = q ⊗ x ⊗ q * , as we explain later in this document).

CAUTION: Not all quaternion definitions are the same. Some authors write the products as ib instead of bi, and therefore they get the property k = ji = -ij, which results in ijk = 1 and a left-handed quaternion. Also, many authors place the real part at the end position, yielding Q = ia+jb+kc+d. These choices have no fundamental implications but make the whole formulation different in the details. Please refer to Section 3 for further explanations and disambiguation.

CAUTION: There are additional conventions that also make the formulation different in details. They concern the "meaning" or "interpretation" we give to the rotation operators, either rotating vectors or rotating reference frames -which, essentially, constitute opposite operations. Refer also to Section 3 for further explanations and disambiguation. NOTE: Among the different conventions exposed above, this document concentrates on the Hamilton convention, whose most remarkable property is the definition (2). A proper and grounded disambiguation requires to first develop a significant amount of material; therefore, this disambiguation is relegated to the aforementioned Section 3.

Alternative representations of the quaternion

The real + imaginary notation {1, i, j, k} is not always convenient for our purposes. Provided that the algebra (2) is used, a quaternion can be posed as a sum scalar + vector,

Q = q w + q x i + q y j + q z k ⇔ Q = q w + q v , (5)
where q w is referred to as the real or scalar part, and q v = q x i + q y j + q z k = (q x , q y , q z ) as the imaginary or vector part. 1 It can be also defined as an ordered pair scalar-vector

Q = q w , q v . (6)
We mostly represent a quaternion Q as a 4-vector q ,

q q w q v =     q w q x q y q z     , (7)
which allows us to use matrix algebra for operations involving quaternions. At certain occasions, we may allow ourselves to mix notations by abusing of the sign "=". Typical examples are real quaternions and pure quaternions, general: q = q w + q v = q w q v ∈ H , real: q w = q w 0 v ∈ R , pure:

q v = 0 q v ∈ H p . (8)
1.2 Main quaternion properties

Sum

The sum is straightforward,

p q = p w p v q w q v = p w q w p v q v . ( 9
)
By construction, the sum is commutative and associative, p + q = q + p (10) p + (q + r) = (p + q) + r .

(11)

Product

Denoted by ⊗, the quaternion product requires using the original form (1) and the quaternion algebra (2). Writing the result in vector form gives

p ⊗ q =    
p w q wp x q xp y q yp z q z p w q x + p x q w + p y q zp z q y p w q yp x q z + p y q w + p z q x p w q z + p x q yp y q x + p z q w     .

(

) 12
This can be posed also in terms of the scalar and vector parts,

p ⊗ q = p w q w -p v q v p w q v + q w p v + p v ×q v , (13)
where the presence of the cross-product reveals that the quaternion product is not commutative in the general case,

p ⊗ q = q ⊗ p . ( 14
)
Exceptions to this general non-commutativity are limited to the cases where p v ×q v = 0, which happens whenever one quaternion is real, p = p w or q = q w , or when both vector parts are parallel, p v q v . Only in these cases the quaternion product is commutative.

The quaternion product is however associative,

(p ⊗ q) ⊗ r = p ⊗ (q ⊗ r) , (15)
and distributive over the sum, p ⊗ (q + r) = p ⊗ q + p ⊗ r and (p + q) ⊗ r = p ⊗ r + q ⊗ r .

The product of two quaternions is bi-linear and can be expressed as two equivalent matrix products, namely

q 1 ⊗ q 2 = [q 1 ] L q 2 and q 1 ⊗ q 2 = [q 2 ] R q 1 , (17)
where [q] L and [q] R are respectively the left-and right-quaternion-product matrices, which are derived from ( 12) and ( 17) by simple inspection,

[q] L =     q w -q x -q y -q z q x q w -q z q y q y q z q w -q x q z -q y q x q w     ,

[q] R =    
q w -q x -q y -q z q x q w q z -q y q y -q z q w q x q z q y -q x q w     , (18) or more concisely, from ( 13) and ( 17),

[q] L = q w I + 0 -q v q v [q v ] × , [q] R = q w I + 0 -q v q v -[q v ] × . (19)
Here, the skew operator2 [•] × produces the cross-product matrix,

[a] ×   0 -a z a y a z 0 -a x -a y a x 0   , (20)
which is a skew-symmetric matrix, [a] × = -[a] × , equivalent to the cross product, i.e.,

[a] × b = a×b , ∀ a, b ∈ R 3 . (21)
Finally, since

(q ⊗ x) ⊗ p = [p] R [q] L x and q ⊗ (x ⊗ p) = [q] L [p] R x (22)
we have the relation

[p] R [q] L = [q] L [p] R , (23)
that is, left and right quaternion product matrices commute. Further properties of these matrices are provided in Section 2.8. Quaternions endowed with the product operation ⊗ form a non-commutative group. The group's elements identity, q 1 = 1, and inverse, q -1 , are explored below.

Identity

The identity quaternion q 1 with respect to the product is such that q 1 ⊗ q = q ⊗ q 1 = q. It corresponds to the real product identity '1' expressed as a quaternion,

q 1 = 1 = 1 0 v .

Conjugate

The conjugate of a quaternion is defined by

q * q w -q v = q w -q v . ( 24
)
This has the properties q ⊗ q * = q * ⊗ q = q 2 w + q 2 x + q 2 y + q 2 z =

q 2 w + q 2 x + q 2 y + q 2 z 0 v , (25) and
(p ⊗ q) * = q * ⊗ p * .

(26)

Norm

The norm of a quaternion is defined by q √ q ⊗ q * = √ q * ⊗ q = q 2 w + q 2 x + q 2 y + q 2 z ∈ R .

It has the property p ⊗ q = q ⊗ p = p q .

(28)

Inverse

The inverse quaternion q -1 is such that the quaternion times its inverse gives the identity,

q ⊗ q -1 = q -1 ⊗ q = q 1 . (29)
It can be computed with q -1 = q * / q 2 . (30)

Unit or normalized quaternion

For unit quaternions, q = 1, and therefore

q -1 = q * . ( 31
)
When interpreting the unit quaternion as an orientation specification, or as a rotation operator, this property implies that the inverse rotation can be accomplished with the conjugate quaternion. Unit quaternions can always be written in the form,

q = cos θ u sin θ , (32)
where u = u x i + u y j + u z k is a unit vector and θ is a scalar. From (28), unit quaternions endowed with the product operation ⊗ form a non commutative group, where the inverse coincides with the conjugate.

Additional quaternion properties 1.3.1 Quaternion commutator

The quaternion commutator is defined as [p, q] p ⊗ qq ⊗ p. We have from (13), p ⊗ qq ⊗ p = 2 p v ×q v .

(33)

This has as a trivial consequence,

p v ⊗ q v -q v ⊗ p v = 2 p v ×q v . (34)
We will use this property later on.

Product of pure quaternions

Pure quaternions are those with null real or scalar part, Q = q v or q = [0, q v ]. We have from (13),

p v ⊗ q v = -p v q v + p v ×q v = -p v q v p v ×q v . ( 35
)
This implies

q v ⊗ q v = -q v q v = -q v 2 , (36)
and for pure unitary quaternions u ∈ H p , u = 1,

u ⊗ u = -1 , (37)
which is analogous to the standard imaginary case, i • i = -1.

Natural powers of pure quaternions

Let us define q n , n ∈ N, as the n-th power of q using the quaternion product ⊗. Then, if v is a pure quaternion and we let v = u θ, with θ = v ∈ R and u unitary, we get from (36) the cyclic pattern

v 2 = -θ 2 , v 3 = -u θ 3 , v 4 = θ 4 , v 5 = u θ 5 , v 6 = -θ 6 , • • • (38)

Exponential of pure quaternions

The quaternion exponential is a function on quaternions analogous to the ordinary exponential function. Exactly as in the real exponential case, it is defined as the absolutely convergent power series,

e q ∞ k=0 1 k! q k ∈ H . (39)
Clearly, the exponential of a real quaternion coincides exactly with the ordinary exponential function.

More interestingly, the exponential of a pure quaternion v = v x i + v y j + v z k is a new quaternion defined by,

e v = ∞ k=0 1 k! v k ∈ H . ( 40
)
Letting v = u θ, with θ = v ∈ R and u unitary, and considering (38), we group the scalar and vector terms in the series,

e uθ = 1 - θ 2 2! + θ 4 4! + • • • + uθ - uθ 3 3! + uθ 5 5! + • • • (41)
and recognize in them, respectively, the series of cos θ and sin θ. 3 This results in

e v = e u θ = cos θ + u sin θ = cos θ u sin θ , (42)
which constitutes a beautiful extension of the Euler formula, e iθ = cos θ + i sin θ, defined for imaginary numbers. Notice that since e v 2 = cos 2 θ + sin 2 θ = 1, the exponential of a pure quaternion is a unit quaternion. Notice also the property,

e -v = (e v ) * . ( 43
)
For small angle quaternions we avoid the division by zero in u = v/ v by expressing the Taylor series of sin θ and cos θ and truncating, obtaining varying degrees of the approximation,

e v ≈ 1 -θ 2 /2 v 1 -θ 2 /6 ≈ 1 v --→ θ→0 1 0 . ( 44
)

Exponential of general quaternions

Due to the non-commutativity property of the quaternion product, we cannot write for general quaternions p and q that e p+q = e p e q . However, commutativity holds when any of the product members is a scalar, and therefore, e q = e qw+qv = e qw e qv .

(45)

Then, using (42) with uθ = q v we get e q = e qw cos q v qv qv sin q v . (46)

1.3.6 Logarithm of unit quaternions It is immediate to see that, if q = 1, log q = log(cos θ + u sin θ) = log(e u θ ) = u θ = 0 u θ , (47)
that is, the logarithm of a unit quaternion is a pure quaternion. The angle-axis values are obtained easily by inverting (42), u = q v / q v (48) θ = arctan( q v , q w ) .

(49)

For small angle quaternions, we avoid division by zero by expressing the Taylor series of arctan(x) and truncating,4 obtaining varying degrees of the approximation,

log(q) = uθ = q v arctan( q v , q w ) q v ≈ q v q w 1 - q v 2 3q 2 w ≈ q v --→ θ→0 0 . (50)

Logarithm of general quaternions

By extension, if q is a general quaternion, log q = log( q q q ) = log q + log q q = log q + u θ = log q u θ .

(51)

1.3.8 Exponential forms of the type q t

We have, for q ∈ H and t ∈ R, q t = exp(log(q t )) = exp(t log(q)) .

If q = 1, we can write q = [cos θ, u sin θ], thus log(q) = uθ, which gives

q t = exp(t uθ) = cos tθ u sin tθ . ( 53
)
Because the exponent t has ended up as a linear multiplier of the angle θ, it can be seen as a linear angular interpolator. We will develop this idea in Section 2.7.

2 Rotations and cross-relations 2.1 The 3D vector rotation formula We illustrate in Fig. 1 the rotation, following the right-hand rule, of a general 3D vector x, by an angle φ, around the axis defined by the unit vector u. This is accomplished by decomposing the vector x into a part x || parallel to u, and a part x ⊥ orthogonal to u, so that

x || x ? x 0 x 0 ? u x x ? = x u u > x x || = u u > x x = x || + x ? ↵ Figure
x = x || + x ⊥ .
These parts can be computed easily (α is the angle between the vector x and the axis u),

x || = u ( x cos α) = u u x x ⊥ = x -x || = x -u u x .
Upon rotation, the parallel part does not rotate,

x || = x || ,
and the orthogonal part experiences a planar rotation in the plane normal to u. That is, if we create an orthogonal base {e 1 , e 2 } of this plane with

e 1 = x ⊥ e 2 = u×x ⊥ = u×x ,
satisfying e 1 = e 2 , then x ⊥ = e 1 •1 + e 2 •0. A rotation of φ rad on this plane produces,

x ⊥ = e 1 cos φ + e 2 sin φ , which develops as,

x ⊥ = x ⊥ cos φ + (u×x) sin φ .
Adding the parallel part yields the expression of the rotated vector, x = x || + x ⊥ , which is known as the vector rotation formula,

x = x || + x ⊥ cos φ + (u × x) sin φ .
(54)

The rotation group SO(3)

In R 3 , the rotation group SO(3) is the group of rotations around the origin under the operation of composition. Rotations are linear transformations that preserve vector length and relative vector orientation (i.e., handedness). Its importance in robotics is that it represents rotations of rigid bodies in 3D space: a rigid motion requires precisely that distances, angles and relative orientations within a rigid body be preserved upon motion -otherwise, if norms, angles or relative orientations are not kept, the body could not be considered rigid.

Let us then define rotations through an operator that satisfies these properties. A rotation operator r : R 3 → R 3 ; v → r(v) acting on vectors v ∈ R 3 can be defined from the metrics of Euclidean space, constituted by the dot and cross products, as follows.

• Rotation preserves the vector norm,

r(v) = r(v), r(v) = v, v v , ∀v ∈ R 3 . (55a)
• Rotation preserves angles between vectors,

r(v), r(w) = v, w = v w cos α , ∀v, w ∈ R 3 . ( 55b
)
• Rotation preserves the relative orientations of vectors,

u × v = w ⇐⇒ r(u) × r(v) = r(w) . ( 56
)
It is easily proved that the first two conditions are equivalent. We can thus define the rotation group SO(3) as,

SO(3) : {r : R 3 → R 3 / ∀ v, w ∈ R 3 , r(v) = v , r(v)×r(w) = r(v×w)} . ( 57
)
The rotation group is typically represented by the set of rotation matrices. However, quaternions constitute also a good representation of it. The aim of this chapter is to

Rotation matrix, R

Quaternion, q

Parameters 3 × 3 = 9 1 + 3 = 4 Degrees of freedom 3 3 Constraints 9 -3 = 6 4 -3 = 1 Constraints RR = I ; det(R) = +1 q ⊗ q * = 1 ODE Ṙ = R [ω] × q = 1 2 q ⊗ ω Exponential map R = exp([uφ] × ) q = exp(uφ/2) Logarithmic map log(R) = [uφ] × log(q) = uφ/2 Relation to SO(3) Single cover Double cover Identity I 1 Inverse R q * Composition R 1 R 2 q 1 ⊗ q 2 Rotation operator R = I + sin φ [u] × + (1 -cos φ) [u] 2 × q = cos φ/2 + u sin φ/2 Rotation action R x q ⊗ x ⊗ q * Interpolation R t = I + sin tφ [u] × + (1-cos tφ) [u] 2 × q t = cos tφ/2 + u sin tφ/2 R 1 (R 1 R 2 ) t q 1 ⊗ (q * 1 ⊗ q 2 ) t q 1 sin((1-t)∆θ) sin(∆θ) + q 2 sin(t∆θ) sin(∆θ)
Cross relations

R{q} = (q 2 w -q v q v ) I + 2 q v q v + 2 q w [q v ] × R{-q} = R{q} double cover R{1} = I identity R{q * } = R{q} inverse R{q 1 ⊗ q 2 } = R{q 1 } R{q 2 } composition R{q t } = R{q} t interpolation
show that both representations are equally valid. They exhibit a lot of similarities, both conceptually and algebraically, as the reader will appreciate in Table 1. Perhaps, the most important difference is that the unit quaternion group constitutes a double cover of SO(3) (thus technically not being SO(3) itself), something that is not critical in most of our applications. 5 The table is inserted upfront for the sake of a rapid comparison and evaluation. The rotation matrix and quaternion representations of SO(3) are explored in the following sections.

The rotation group and the rotation matrix

The operator r() is linear, since it is defined from the scalar and vector products, which are linear. It can therefore be represented by a matrix R ∈ R 3×3 , which produces rotations to vectors v ∈ R 3 through the matrix product,

r(v) = R v . ( 58
)
Injecting it in (55a), using the dot product a, b = a b and developing we have that for all v,

(Rv) (Rv) = v R Rv = v v , (59)
yielding the orthogonality condition on R,

R R = I = R R . ( 60
)
The condition above is indeed a condition of orthogonality, since we can observe from it that, by writing R = [r 1 , r 2 , r 3 ] and substituting above, the column vectors r i of R, with i ∈ {1, 2, 3}, are of unit length and orthogonal to each other,

r i , r i = r i r i = 1 r i , r j = r i r j = 0 , if i = j .
The set of transformations keeping vector norms and angles is for this reason called the Orthogonal group, denoted O(3). The orthogonal group includes rotations (which are rigid motions) and reflections (which are not rigid). The notion of group here means essentially (and informally) that the product of two orthogonal matrices is always an orthogonal matrix,6 and that each orthogonal matrix admits an inverse. In effect, the orthogonality condition (60) implies that the inverse rotation is achieved with the transposed matrix,

R -1 = R . ( 61
)
Adding the relative orientations condition (56) guarantees rigid body motion (hence discarding reflections), and results in one additional constraint on R,7 det(R) = 1 .

(62)

Orthogonal matrices with positive unit determinant are commonly referred to as proper or special. The set of such special orthogonal matrices is a subgroup of O(3) named the Special Orthogonal group SO(3). Being a group, the product of two rotation matrices is always a rotation matrix.8

The exponential map

The exponential map (and the logarithmic map, which we see in the next section) is a powerful mathematical tool for working in the rotational 3D space with ease and rigor. It represents the entrance door to a corpus of infinitesimal calculus suited for the rotational space. The exponential map allows us to properly define derivatives, perturbations, and velocities, and to manipulate them. It is therefore essential in estimation problems in the space of rotations or orientations. Rotations constitute rigid motions. This rigidity implies that it is possible to define a continuous trajectory or path r(t) in SO(3) that continuously rotates the rigid body from its initial orientation, r(0), to its current orientation, r(t). Being continuous, it is legitimate to investigate the time-derivatives of such transformations. We do so by deriving the properties (60) and (62) that we have just seen.

First of all, we notice that it is impossible to continuously escape the unit determinant condition (62) while satisfying (60), because this would imply a jump of the determinant from +1 to -1.9 Therefore we only need to investigate the time-derivative of the orthogonality condition (60). This reads

d dt (R R) = Ṙ R + R Ṙ = 0 , (63)
which results in

R Ṙ = -(R Ṙ) , (64)
meaning that the matrix R Ṙ is skew-symmetric (i.e., it is equal to the negative of its transpose). The set of skew-symmetric 3 × 3 matrices is denoted so(3), and receives the name of the Lie algebra of SO(3). Skew-symmetric 3 × 3 matrices have the form, they have 3 DOF, and correspond to cross-product matrices, as we introduced already in (20). This establishes a one-to-one mapping ω ∈ R 3 ↔ [ω] × ∈ so(3). Let us then take a vector ω = (ω x , ω y , ω z ) ∈ R 3 and write

[ω] ×   0 -ω z ω y ω z 0 -ω x -ω y ω x 0   ; (65) v 2 R 3 [v] ⇥ 2 so(3) R 2 SO(3) Exp(•) exp(•) [•] ⇥
R Ṙ = [ω] × . ( 66
)
This leads to the ordinary differential equation (ODE),

Ṙ = R [ω] × . (67)
Around the origin, we have R = I and the equation above reduces to Ṙ = [ω] × . Thus, we can interpret the Lie algebra so(3) as the space of the derivatives of r(t) at the origin; it constitutes the tangent space to SO(3), or the velocity space. Following these facts, we can very well call ω the vector of instantaneous angular velocities.

If ω is constant, the differential equation above can be time-integrated as

R(t) = R(0) e [ω] × t = R(0) e [ωt] × (68)
where the exponential e [x] × is defined by its Taylor series, as we see in the following section. Since R(0) and R(t) are rotation matrices, then clearly e [ωt] × = R(0) R(t) is a rotation matrix. Defining the vector v ω∆t as the rotation vector encoding the full rotation over the period ∆t, we have

R = e [v] × . ( 69
)
This is known as the exponential map, an application from so(3) to SO(3),

exp : so(3) → SO(3) ; [v] × → exp([v] × ) = e [v] × . ( 70
)

The capitalized exponential map

The exponential map above is sometimes expressed with some abuse of notation, i.e., confounding v ∈ R 3 with [v] × ∈ so(3). To avoid possible ambiguities, we opt for writing this new application R 3 → SO(3) with an explicit notation using a capitalized Exp, having (see Fig. 2)

Exp : R 3 → SO(3) ; v → Exp(v) = e [v] × . ( 71
)
Its relation with the exponential map is trivial,

Exp(v) exp([v] × ) . (72)
In the following sections we'll see that the vector v, called the rotation vector or the angle-axis vector, encodes through v = ω∆t = φu the angle φ and axis u of rotation.

Rotation matrix and rotation vector: the Rodrigues rotation formula

The rotation matrix is defined from the rotation vector v = φu through the exponential map (69), with the cross-product matrix [v] × = φ [u] × as defined in (20). The Taylor expansion of (69) with v = φu reads,

R = e φ[u] × = I + φ [u] × + 1 2 φ 2 [u] 2 × + 1 3! φ 3 [u] 3 × + 1 4! φ 4 [u] 4 × + . . . (73)
When applied to unit vectors, u, the matrix

[u] × satisfies [u] 2 × = uu -I (74) [u] 3 × = -[u] × , (75)
and thus all powers of [u] × can be expressed in terms of

[u] × and [u] 2 × in a cyclic pattern, [u] 4 × = -[u] 2 × [u] 5 × = [u] × [u] 6 × = [u] 2 × [u] 7 × = -[u] × • • • . ( 76
)
Then, grouping the Taylor series in terms of [u] × and [u] 2 × , and identifying in them, respectively, the series of sin φ and cos φ, leads to a closed form to obtain the rotation matrix from the rotation vector, the so called Rodrigues rotation formula,

R = I + sin φ [u] × + (1 -cos φ) [u] 2 × , (77)
which we denote R{v} Exp(v). This formula admits some variants, e.g., using (74),

R = I cos φ + [u] × sin φ + uu (1 -cos φ) . (78)

The logarithmic maps

We define the logarithmic map as the inverse of the exponential map, log :

SO(3) → so(3) ; R → log(R) = [u φ] × , (79)
with

φ = arccos trace(R) -1 2 (80) u = (R -R ) ∨ 2 sin φ , (81)
where

• ∨ is the inverse of [•] × , that is, ([v] × ) ∨ = v and [V ∨ ] × = V.
We also define a capitalized version Log, which allows us to recover the rotation vector v = uφ ∈ R 3 directly from the rotation matrix,

Log : SO(3) → R 3 ; R → Log(R) = u φ . ( 82a
)
Its relation with the logarithmic map is trivial,

Log(R) (log(R)) ∨ . ( 83
)

The rotation action

Rotating a vector x by an angle φ around the unit axis u is performed with the linear product

x = R x , (84)
where R = Exp(uφ). This can be shown by developing (84), using ( 77), ( 74) and ( 75),

x = R x = (I + sin φ [u] × + (1 -cos φ) [u] 2 × ) x = x + sin φ [u] × x + (1 -cos φ) [u] 2 × x = x + sin φ(u×x) + (1 -cos φ)(uu -I) x = x + x ⊥ + sin φ(u×x) -(1 -cos φ) x ⊥ = x + (u×x) sin φ + x ⊥ cos φ , (85)
which is precisely the vector rotation formula (54).

The rotation group and the quaternion

For didactical purposes, we are interested in highlighting the connections between quaternions and rotation matrices as representations of the rotation group SO(3). For this, the well-known formula of the quaternion rotation action, which reads,

r(v) = q ⊗ v ⊗ q * , (86)
is here taken initially as an hypothesis. This allows us to develop the full quaternion section with a discourse that retraces the one we used for the rotation matrix. The exactness of this hypothesis will be proved a little later, in Section 2.4.5, thus validating the approach.

Let us then inject the rotation above into the orthogonality condition (55a), and develop it using (28) as

q ⊗ v ⊗ q * = q 2 v = v . ( 87
)
This yields q 2 = 1, that is, the unit norm condition on the quaternion, which reads,

q * ⊗ q = 1 = q ⊗ q * . ( 88
)
This condition is akin to the one we encountered for rotation matrices, see (60), which reads R R = I = RR . We encourage the reader to stop at their similarities for a second.

Similarly, we show that the relative orientation condition ( 56) is satisfied by construction (we use (34) twice),

r(v) × r(w) = (q ⊗ v ⊗ q * ) × (q ⊗ w ⊗ q * ) = 1 2 (q ⊗ v ⊗ q * ) ⊗ (q ⊗ w ⊗ q * ) -(q ⊗ w ⊗ q * ) ⊗ (q ⊗ v ⊗ q * ) = 1 2 (q ⊗ v ⊗ w ⊗ q * -q ⊗ w ⊗ v ⊗ q * ) = 1 2 (q ⊗ (v ⊗ w -w ⊗ v) ⊗ q * ) = q ⊗ (v × w) ⊗ q * = r(v × w) . ( 89
)
The set of unit quaternions forms a group under the operation of multiplication. This group is topologically a 3-sphere, that is, the 3-dimensional surface of the unit sphere of R 4 , and is commonly noted as S 3 .

The exponential map

Let us consider a unit quaternion q ∈ S 3 , that is, q * ⊗ q = 1, and let us proceed as we did for the orthogonality condition of the rotation matrix, R R = I. Taking the time derivative,

d(q * ⊗ q) dt = q * ⊗ q + q * ⊗ q = 0 , (90)
it follows that

q * ⊗ q = -( q * ⊗ q) = -(q * ⊗ q) * , (91)
which means that q * ⊗ q is a pure quaternion (i.e., it is equal to minus its conjugate, therefore its real part is zero). We thus take a pure quaternion Ω ∈ H p and write,

q * ⊗ q = Ω = 0 Ω ∈ H p . (92)
Left-multiplication by q yields the differential equation, Around the origin, we have q = 1 and the equation above reduces to q = Ω ∈ H p . Thus, the space H p of pure quaternions constitutes the tangent space, or the Lie Algebra, of the unit sphere S 3 of quaternions. In the quaternion case, however, this space is not directly the velocity space, but rather the space of the half-velocities, as we will see soon.

q = q ⊗ Ω . ( 93
) v 2 R 3 Exp(•) exp(•) ÷2 V 2 H p q 2 S 3
If Ω is constant, the differential equation can be integrated as

q(t) = q(0) ⊗ e Ω t , (94)
where, since q(0) and q(t) are unit quaternions, the exponential e Ωt is also a unit quaternion -something we already knew from the quaternion exponential (42). Defining V Ω∆t we have

q = e V . (95)
This is again an exponential map: an application from the space of pure quaternions to the space of rotations represented by unit quaternions, exp :

H p → S 3 ; V → exp(V) = e V (96)

The capitalized exponential map

As we will see, the pure quaternion V in the exponential map (96) encodes, through V = θu = φu/2, the axis of rotation u and the half of the rotated angle, θ = φ/2. We will provide ample explanations to this half-angle fact very soon, mainly in Sections 2.4.5, 2.4.6 and 2.8. By now, let it suffice to say that, since the rotation action is accomplished by the double product x = q ⊗ x ⊗ q * , the vector x experiences a rotation which is 'twice' the one encoded in q, or equivalently, the quaternion q encodes 'half' the intended rotation on x.

In order to express a direct relation between the angle-axis rotation parameters, v = φu ∈ R 3 , and the quaternion, we define a capitalized version of the exponential map, which captures the half-angle effect (see Fig. 3),

Exp : R 3 → S 3 ; v → Exp(v) = e v/2
(97)

Its relation to the exponential map is trivial,

Exp(v) exp(v/2) . ( 98
)
It is also convenient to introduce the vector of angular velocities ω = 2Ω ∈ R 3 , so that (93) and (94) become, q = 1 2 q ⊗ ω (99)

q = e ωt/2 . (100)

Quaternion and rotation vector

Let v = φu be a rotation vector representing a rotation of φ rad around the axis u. Then, the exponential map can be developed using an extension of the Euler formula (see (38-42) for a complete development),

q Exp(φu) = e φu/2 = cos φ 2 + u sin φ 2 = cos(φ/2) u sin(φ/2) . ( 101
)
We call this the rotation vector to quaternion conversion formula, and will be denoted in this document by q = q{v} Exp(v).

The logarithmic maps

We define the logarithmic map as the inverse of the exponential map, log :

S 3 → H p ; q → log(q) = uθ , (102)
which is of course the definition we gave for the quaternion logarithm in Section 1.3.6. We also define the capitalized logarithmic map, which directly provides the angle φ and axis u of rotation in Cartesian 3-space, Log : S 3 → R 3 ; q → Log(q) = uφ .

(103)

Its relation with the logarithmic map is trivial,

Log(q) 2 log(q) . ( 104
)
For its implementation we use the 4-quadrant version of arctan(y, x). From (101),

φ = 2 arctan( q v , q w ) (105a) u = q v / q v . ( 105b
)
For small-angle quaternions, (105b) diverges. We then use the a truncated Taylor series for the arctan() function, getting,

Log(q) = θu ≈ 2 q v q w 1 - q v 2 3q 2 w . (106)

The rotation action

We are finally in the position of proving our hypothesis (86) for the vector rotation using quaternions, thus validating all the material presented so far. Rotating a vector x by an angle φ around the axis u is performed with the double quaternion product, also known as the sandwich product,

x = q ⊗ x ⊗ q * , (107)
where q = Exp(uφ), and where the vector x has been written in quaternion form, that is,

x = xi + yj + zk = 0 x ∈ H p . ( 108
)
To show that this double product does perform the desired vector rotation, we use ( 13), (101), and basic vector and trigonometric identities, to develop (107) as follows,

x = q ⊗ x ⊗ q * = cos φ 2 + u sin φ 2 ⊗ (0 + x) ⊗ cos φ 2 -u sin φ 2 = x cos 2 φ 2 + (u ⊗ x -x ⊗ u) sin φ 2 cos φ 2 -u ⊗ x ⊗ u sin 2 φ 2 = x cos 2 φ 2 + 2(u×x) sin φ 2 cos φ 2 -(x(u u) -2u(u x)) sin 2 φ 2 = x(cos 2 φ 2 -sin 2 φ 2 ) + (u×x)(2 sin φ 2 cos φ 2 ) + u(u x)(2 sin 2 φ 2 ) = x cos φ + (u×x) sin φ + u(u x)(1 -cos φ) = (x -u u x) cos φ + (u×x) sin φ + u u x = x ⊥ cos φ + (u×x) sin φ + x || , (109)
which is precisely the vector rotation formula (54).

2.4.6

The double cover of the manifold of SO(3).

Consider a unit quaternion q. When regarded as a regular 4-vector, the angle θ between q and the identity quaternion q 1 = [1, 0, 0, 0] representing the origin of orientations is,

cos θ = q 1 q = q(1) = q w . ( 110
)
At the same time, the angle φ rotated by the quaternion q on objects in 3D space satisfies

q = q w q v = cos φ/2 u sin φ/2 . ( 111
)
That is, we have q w = cos θ = cos φ/2, so the angle between a quaternion vector and the identity in 4D space is half the angle rotated by the quaternion in 3D space,

θ = φ/2 . ( 112
)
We illustrate this double cover in Fig. 4. By the time the angle between the two quaternion vectors is θ = π/2, the 3D rotation has already achieved φ = π, which is half a turn. And by the time the quaternion vector has made a half turn, θ = π, the 3D rotation has completed a full turn. The second half turn of the quaternion vector, π < θ < 2π, represents a second full turn of the 3D rotation, 2π < φ < 4π, that is, a second cover of the rotation manifold.

✓ q ⇡ ✓ q 1 x 0 = 2✓ u x
Figure 4: Double cover of the rotation manifold. Left: the quaternion q in the unit 3sphere defines an angle θ with the identity quaternion q 1 . Center: the resulting 3D rotation x = q ⊗ x ⊗ q * has double angle φ than that of the original quaternion. Right: Superposing the 4D and 3D rotation planes, observe how one turn of the quaternion q over the 3-sphere (red) represents two turns of the rotated vector x in 3D space (blue).

Rotation matrix and quaternion

As we have just seen, given a rotation vector v = u φ, the exponential maps for the unit quaternion and the rotation matrix produce rotation operators q = Exp(u φ) and R = Exp(u φ) that rotate vectors x exactly the same angle φ around the same axis u.10 That is, if

∀v, x ∈ R 3 , q = Exp(v), R = Exp(v) (113) then, q ⊗ x ⊗ q * = R x . ( 114
)
As both sides of this identity are linear in x, an expression of the rotation matrix equivalent to the quaternion is found by developing the left hand side and identifying terms on the right, yielding the quaternion to rotation matrix formula,

R =   q 2 w + q 2
xq 2 yq 2 z 2(q x q yq w q z ) 2(q x q z + q w q y ) 2(q x q y + q w q z ) q 2 wq 2 x + q 2 yq 2 z 2(q y q zq w q x ) 2(q x q zq w q y ) 2(q y q z + q w q x ) q 2 w -

q 2 x -q 2 y + q 2 z   , (115)
denoted throughout this document by R = R{q}. The matrix form of the quaternion product (17-19) provides us with an alternative formula, since

q ⊗ x ⊗ q * = [q * ] R [q] L 0 x = 0 R x , (116)
which leads after some easy developments to

R = (q 2 w -q v q v ) I + 2 q v q v + 2 q w [q v ] × . (117)
The rotation matrix R has the following properties with respect to the quaternion,

R{[1, 0, 0, 0] } = I (118) R{-q} = R{q} (119) R{q * } = R{q} (120) R{q 1 ⊗ q 2 } = R{q 1 }R{q 2 } , (121)
where we observe that: (118) the identity quaternion encodes the null rotation; (119) a quaternion and its negative encode the same rotation, defining a double cover of SO(3);

(120) the conjugate quaternion encodes the inverse rotation; and (121) the quaternion product composes consecutive rotations in the same order as rotation matrices do. Additionally, we have the property

R{q t } = R{q} t , (122)
which relates the spherical interpolations of the quaternion and rotation matrix over a running scalar t.

Rotation composition

Quaternion composition is done similarly to rotation matrices, i.e., with appropriate quaternion-and matrix-products, and in the same order (Fig. 5),

q AC = q AB ⊗ q BC , R AC = R AB R BC . ( 123
)
This comes immediately from the associative property of the involved products,

x A = q AB ⊗ x B ⊗ q * AB x A = R AB x B = q AB ⊗ (q BC ⊗ x C ⊗ q * BC ) ⊗ q * AB = R AB (R BC x C ) = (q AB ⊗ q BC ) ⊗ x C ⊗ (q * BC ⊗ q * AB ) = (R AB R BC ) x C = (q AB ⊗ q BC ) ⊗ x C ⊗ (q AB ⊗ q BC ) * = R AC x C . = q AC ⊗ x C ⊗ q * AC , A B C q AC q AB q BC q AC = q AB ⌦ q BC Figure 5: Rotation composition.
In R 2 , we would simply do θ AC = θ AB + θ BC , with an operation 'sum' that is commutative. In R 3 composition satisfies q AC = q AB ⊗ q BC and, in matrix form, R AC = R AB R BC . These operators are not commutative and one must respect the order strictly -a proper notation helps: 'AB' chains with 'BC' to create 'AC'.

A comment on notation A proper notation helps determining the right order of the factors in the composition, especially for compositions of several rotations (see Fig. 5). For example, let q ji (resp. R ji ) represent a rotation from situation i to situation j, that is,

x j = q ji ⊗ x i ⊗ q * ji (resp. x j = R ji x i ).
Then, given a number of rotations represented by the quaternions q OA , q AB , q BC , q OX , q XZ , we just have to chain the indices and get:

q OC = q OA ⊗ q AB ⊗ q BC R OC = R OA R AB R BC ,
and knowing that the opposite rotation corresponds to the conjugate, q ji = q * ij , or the transpose, R ji = R ij , we also have

q ZA = q * XZ ⊗ q * OX ⊗ q OA R ZA = R XZ R OX R OA = q ZX ⊗ q XO ⊗ q OA = R ZX R XO R OA .

Spherical linear interpolation (SLERP)

Quaternions are very handy for computing proper orientation interpolations. Given two orientations represented by quaternions q 0 and q 1 , we want to find a quaternion function q(t), t ∈ [0, 1], that linearly interpolates from q(0) = q 0 to q(1) = q 1 . This interpolation is such that, as t evolves from 0 to 1, a body will continuously rotate from orientation q 0 to orientation q 1 , at constant speed along a fixed axis.

Method 1 A first approach uses quaternion algebra, and follows a geometric reasoning in R 3 that should be easily related to the material presented so far. First, compute the orientation increment ∆q from q 0 to q 1 such that q 1 = q 0 ⊗ ∆q,

∆q = q * 0 ⊗ q 1 . ( 124
)
Then obtain the associated rotation vector, ∆φ = u∆φ, using the logarithmic map,11

u ∆φ = Log(∆q) . (125)
Finally, keep the rotation axis u and take a linear fraction of the rotation angle, δφ = t∆φ.

Put it in quaternion form through the exponential map, δq = Exp(u δφ), and compose it with the original quaternion to get the interpolated result,

q(t) = q 0 ⊗ Exp(t u ∆φ) . ( 126
)
The whole process can be written as q(t) = q 0 ⊗ Exp(t Log(q * 0 ⊗ q 1 )), which reduces to

q(t) = q 0 ⊗ (q * 0 ⊗ q 1 ) t , (127)
and which is usually implemented (see ( 53)) as,

q(t) = q 0 ⊗ cos(t ∆φ/2) u sin(t ∆φ/2) . ( 128
)
Note: An analogous procedure may be used to define Slerp for rotation matrices, yielding

R(t) = R 0 Exp(t Log(R 0 R 1 )) = R 0 (R 0 R 1 ) t , (129)
where the matrix exponential R t can be implemented using Rodrigues (77), leading to Method 2 Other approaches to Slerp can be developed that are independent of the inners of quaternion algebra, and even independent of the dimension of the space in which the arc is embedded. In particular, see Fig. 6, we can treat quaternions q 0 and q 1 as two unit vectors in the unit sphere, and interpolate in this same space. The interpolated q(t) is the unit vector that follows at a constant angular speed the shortest spherical path joining q 0 to q 1 . This path is the planar arc resulting from intersecting the unit sphere with the plane defined by q 0 , q 1 and the origin (dashed circumference in the figure). For a proof that these approaches are equivalent to the above, see Dam et al. (1998).

R(t) = R 0 I + sin(t∆φ) [u] × + (1 -cos(t∆φ)) [u] 2 × . ( 130
) q 0 q 1 q ? q(t) t ✓ q 0 q 1 q(t) ⇡ ⇡ ✓
The first of these approaches uses vector algebra and follows literally the ideas above. Consider q 0 and q 1 as two unit vectors; the angle12 between them is derived from the scalar product, cos(∆θ) = q 0 q 1 ∆θ = arccos(q 0 q 1 ) .

We proceed as follows. We identify the plane of rotation, that we name here π, and build its ortho-normal basis {q 0 , q ⊥ }, where q ⊥ comes from ortho-normalizing q 1 against q 0 , q ⊥ = q 1 -(q 0 q 1 )q 0 q 1 -(q 0 q 1 )q 0 , (132) so that (see Fig. 6 -right)

q 1 = q 0 cos ∆θ + q ⊥ sin ∆θ . (133)
Then, we just need to rotate q 0 a fraction of the angle, t∆θ, over the plane π, to yield the spherical interpolation, q(t) = q 0 cos(t∆θ) + q ⊥ sin(t∆θ) .

(134)

Method 3 A similar approach, credited to Glenn Davis in Shoemake (1985), draws from the fact that any point on the great arc joining q 0 to q 1 must be a linear combination of its ends (since the three vectors are coplanar). Having computed the angle ∆θ using (131), we can isolate q ⊥ from ( 133) and inject it in (134). Applying the identity sin(∆θ -t∆θ) = sin ∆θ cos t∆θcos ∆θ sin t∆θ, we obtain the Davis' formula (see Eberly (2010) for an alternative derivation),

q(t) = q 0 sin((1 -t)∆θ) sin(∆θ) + q 1 sin(t∆θ) sin(∆θ) . ( 135
)
This formula has the benefit of being symmetric: defining the reverse interpolator s = 1t yields

q(s) = q 1 sin((1 -s)∆θ) sin(∆θ) + q 0 sin(s∆θ) sin(∆θ) .
which is exactly the same formula with the roles of q 0 and q 1 swapped. All these quaternion-based SLERP methods require some care to ensure proper interpolation along the shortest path, that is, with rotation angles φ ≤ π. Due to the quaternion double cover of SO(3) (see Section 2.4.6) only the interpolation between quaternions in acute angles ∆θ ≤ π/2 is done following the shortest path (Fig. 7). Testing for this situation and solving it is simple: if cos(∆θ) = q 0 q 1 < 0, then replace e.g. q 1 by -q 1 and start over.

Quaternion and isoclinic rotations: explaining the magic

This section provides geometrical insights to the two intriguing questions about quaternions, what we call the 'magic':

• How is it that the product q ⊗ x ⊗ q * rotates the vector x?

• Why do we need to consider half-angles when constructing the quaternion through q = e v/2 = [cos φ/2, u sin φ/2]?

We want a geometrical explanation, that is, some rationale that goes beyond the algebraic demonstration (109) and the double cover facts in Section 2.4.6.

To start, let us reproduce here equation ( 116) expressing the quaternion rotation action through the quaternion product matrices [q] L and [q * ] R , defined in (19),

q ⊗ x ⊗ q * = [q * ] R [q] L 0 x = 0 R x . q 0 q 1 q 0 (t) ⇡ q 0 1 = q 1 ✓ > ⇡/2 ✓ 0 < ⇡/2 x 0 x 1 x(t) x 0 (t) > ⇡ 0 < ⇡ q(t)
Figure 7: Ensuring Slerp along the shortest path between the orientations represented by q 0 and q 1 . Left: quaternion rotation plane in 4D space, showing initial and final orientation quaternions, and two possible interpolations, q(t) from q 0 to q 1 , and q (t) from q 0 to -q 1 . Right: vector rotation plane in 3D space: since q 1 = -q 1 , we have

x 1 = q 1 ⊗ x 0 ⊗ q * 1 = q 1 ⊗ x 0 ⊗ q *
1 , that is, both quaternions produce the same rotation. However, the interpolated quaternion q(t) produces the vector x(t) which takes the long path from x 0 to x 1 , while the corrected q 1 = -q 1 yields q (t), producing the vector x (t) along the shortest path from x 0 to x 1 .

For unit quaternions q, the quaternion product matrices [q] L and [q * ] R satisfy two remarkable properties,

[q] [q] = I 4 (136) det([q]) = +1 , (137)
and are therefore elements of SO(4), that is, proper rotation matrices in the R 4 space. To be more specific, they represent a particular type of rotation, named isoclinic rotation, as we explain hereafter. Thus, according to (116), a quaternion rotation corresponds to two chained isoclinic rotations in R 4 . In order to explain the insights of quaternion rotation, we need to understand isoclinic rotations in R 4 . For this, we first need to understand general rotations in R 4 . And to understand rotations in R 4 , we need to go back to R 3 , whose rotations are in fact planar rotations. Let us walk all these steps one by one.

Rotations in R 3 : In R 3 , let us consider the rotations of a vector x around an arbitrary axis represented by the vector u -see Fig. 8, and recall Fig. 1. Upon rotation, vectors parallel to the axis of rotation u do not move, and vectors perpendicular to the axis rotate in the plane π perpendicular to the axis. For general vectors x, the two components of the vector in the plane rotate in this plane, while the axial component remains static. . A rotation of a vector x around an axis u describes a circumference in a plane orthogonal to the axis. The component of x parallel to the axis does not move, and is represented by the small red dot on the axis. The sketch on the right illustrates the radically different behaviors of the rotating point on the plane and axis subspaces.

⇡ 1 ⇡ 2 ⇡ 2 ⇡ 1
Figure 9: Rotations in R 4 . Two orthogonal rotations are possible, on two orthogonal planes π 1 and π 2 . Rotations of a vector x (not drawn) in the plane π 1 cause the two components of the vector parallel to this plane (red dot on π 1 ) to describe a circumference (red circle), leaving the other two components in π 2 unchanged (the red dot). Conversely, rotations in the plane π 2 (blue dot on blue circle in π 2 ) leave the components in π 1 unchanged (blue dot). The sketch on the right better illustrates the situation by resigning to draw unrepresentable perspectives in R 4 , which might be misleading.

Rotations in R 4 : In R 4 , see Fig. 9, due to the extra dimension, the one-dimensional axis of rotation in R 3 becomes a new two-dimensional plane. This second plane provides room for a second rotation. Indeed, rotations in R 4 encompass two independent rotations in two orthogonal planes of the 4-space. This means that every 4-vector of each of these planes rotates in its own plane, and that rotations of general 4-vectors with respect to one plane leave unaffected the vector components in the other plane. These planes are for this reason called 'invariant'.

Isoclinic rotations in R 4 : Isoclinic rotations (from Greek, iso: "equal", klinein: "to incline") are those rotations in R 4 where the angles of rotation in the two invariant planes have the same magnitude. Then, when the two angles have also the same sign,13 we speak of left-isoclinic rotations. And when they have opposite signs, we speak of right-isoclinic rotations. A remarkable property of isoclinic rotations, that we had already seen in ( 23), is that left-and right-isoclinic rotations commute,

⇡ 2 ⇡ 1 ⇡ 2 ⇡ 1 ⇡ 2 ⇡ 1 L R ✓/2 ✓/2 ✓/2 ✓/2
[p] R [q] L = [q] L [p] R . (138)
Quaternion rotations in R 4 and R 3 : Given a unit quaternion q = e u θ/2 , representing a rotation in R 3 of an angle θ around the axis u, the matrix [q] L is a left-isoclinic rotation in R 4 corresponding to the left-multiplication by the quaternion q, and [q * ] R is a rightisoclinic rotation corresponding to the right-multiplication by the quaternion q * . The angles of these isoclinic rotations are exactly of magnitude θ/2,14 and the invariant planes are the same. Then, the rotation expression (116), reproduced once again here,

0 x = q ⊗ x ⊗ q * = [q * ] R [q] L 0 x ,
represents two chained isoclinic rotations to the 4-vector (0, x) , one left-and one right-, each by half the desired rotation angle in R 3 . In one of the invariant planes of R 4 (see Fig. 10), the two half angles cancel out, because they have opposite signs. In the other plane, they sum up to yield the total rotation angle θ. If we define from (116) the resulting rotation matrix, R 4 , one easily realizes that (see also (138)),

R 4 [q * ] R [q] L = [q] L [q * ] R = 1 0 0 R , ( 139
)
where R is the rotation matrix in R 3 , which clearly rotates vectors in the R 3 subspace of R 4 , leaving the fourth dimension unchanged. This discourse is somewhat beyond the scope of the present document. It is also incomplete, for it does not provide, beyond the result in ( 139), an intuition or geometrical explanation for why we need to do q ⊗ x ⊗ q * instead of e.g. q ⊗ x ⊗ q. 15 We include it here just as a means for providing yet another way to interpret rotations by quaternions, with the hope that the reader grasps more intuition about its mechanisms. The interested reader is suggested to consult the appropriate literature on isoclinic rotations in R 4 .

3 Quaternion conventions. My choice.

Quaternion flavors

There are several ways to determine the quaternion. They are basically related to four binary choices:

• The order of its elements -real part first or last:

q = q w q v vs. q = q v q w . (140)
• The multiplication formula -definition of the quaternion algebra:

ij = -ji = k vs. ji = -ij = k , (141a)
which correspond to different handedness, respectively:

right-handed vs. left-handed .

This means that, given a rotation axis u, one quaternion q right {u θ} rotates vectors an angle θ around u using the right hand rule, while the other quaternion q lef t {u θ} uses the left hand rule.

• The function of the rotation operator -rotating frames or rotating vectors:

Passive vs.

Active.

(142)

• In the passive case, the direction of the operation -local-to-global or global-to-local:

x global = q ⊗ x local ⊗ q * vs. x local = q ⊗ x global ⊗ q * (143)

Quaternion type

Hamilton JPL

1 Components order (q w , q v ) (q v , q w ) 2 Algebra ij = k ij = -k Handedness Right-handed Left-handed 3 Function Passive Passive 4 Right-to-left products mean Local-to-Global Global-to-Local Default notation, q q q GL q q LG Default operation x G = q ⊗ x L ⊗ q * x L = q ⊗ x G ⊗ q *
This variety of choices leads to 12 different combinations. Historical developments have favored some conventions over others (Chou, 1992; Yazell, 2009). Today, in the available literature, we find many quaternion flavors such as the Hamilton, the STS 16 , the JPL 17 , the ISS 18 , the ESA 19 , the Engineering, the Robotics, and possibly a lot more denominations. Many of these forms might be identical, others not, but this fact is rarely explicitly stated, and many works simply lack a sufficient description of their quaternion with regard to the four choices above.

These differences impact the respective formulas for rotation, composition, etc., in nonobvious ways. The formulas are thus not compatible, and we need to make a clear choice from the very start.

The two most commonly used conventions, which are also the best documented, are Hamilton (the options on the left in (140-143)) and JPL (the options on the right, with the exception of ( 142)). Table 2 shows a summary of their characteristics. JPL is mostly used in the aerospace domain, while Hamilton is more common to other engineering areas such as robotics -though this should not be taken as a rule.

My choice, that has been taken as early as in equation ( 2), is to take the Hamilton convention, which is right-handed and coincides with many software libraries of widespread use in robotics, such as Eigen, ROS, Google Ceres, and with a vast amount of literature on Kalman filtering for attitude estimation using IMUs (Chou, 1992; Kuipers, 1999; Piniés et al., 2007; Roussillon et al., 2011; Martinelli, 2012, and many others).

The JPL convention is possibly less commonly used, at least in the robotics field. It is extensively described in (Trawny and Roumeliotis, 2005), a reference work that has an aim and scope very close to the present one, but that concentrates exclusively in the JPL convention. The JPL quaternion is used in the JPL literature (obviously) and in key papers by Li, Mourikis, Roumeliotis, and colleagues (see e.g. (Li and Mourikis, 2012;Li et al., 16 Space Transportation System, commonly known as NASA's Space Shuttle. 17 Jet Propulsion Laboratory. 18 International Space Station. 19 European Space Agency.

2014

)), which draw from Trawny and Roumeliotis' document. These works are a primary source of inspiration when dealing with visual-inertial odometry and SLAM -which is what we do.

In the rest of this section we analyze these two quaternion conventions with a little more depth.

Order of the quaternion components

Though not the most fundamental, the most salient difference between Hamilton and JPL quaternions is in the order of the components, with the scalar part being either in first (Hamilton) or last (JPL) position. The implications of such change are quite obvious and should not represent a great challenge of interpretation. In fact, some works with the quaternion's real component at the end (e.g., the C++ library Eigen) are still considered as using the Hamilton convention, as long as the other three aspects are maintained.

We have used the subscripts (w, x, y, z) for the quaternion components for increased clarity, instead of the other commonly used (0, 1, 2, 3). When changing the order, q w will always denote the real part, while it is not clear whether q 0 would also do -in some occasions, one might find things such as q = (q 1 , q 2 , q 3 , q 0 ), with q 0 real and last, but in the general case of q = (q 0 , q 1 , q 2 , q 3 ), the real part at the end would be q 3 . 20 When passing from one convention to the other, we must be careful of formulas involving full 4 × 4 or 3 × 4 quaternion-related matrices, for their rows and/or columns need to be swapped. This is not difficult to do, but it might be difficult to detect and therefore prone to error.

Two curiosities about the components' order are:

• With real part first, the quaternion is naturally interpreted as an extended complex number, of the familiar form real+imaginary. Some of us are comfortable with this representation probably because of this.

• With real part last, the quaternion expressed in vector form, q = x, y, z, w ∈ H, has a format absolutely equivalent to the homogeneous vector in the projective 3D space, p = x, y, z, w ∈ P 3 , where in both cases x, y, z are clearly identified with the three Cartesian axes. When dealing with geometric problems in 3D, this makes the algebra for operating on quaternions and homogeneous vectors more uniform, especially (but not only) if the homogeneous vector is constrained to the unit sphere p = 1.

Specification of the quaternion algebra

The Hamilton convention defines ij = k and therefore,

i 2 = j 2 = k 2 = ijk = -1 , ij = -ji = k , jk = -kj = i , ki = -ik = j , ( 144
)
20 See also footnote 1.

whereas the JPL convention defines ji = k and hence its quaternion algebra becomes,

i 2 = j 2 = k 2 = -ijk = -1 , -ij = ji = k , -jk = kj = i , -ki = ik = j . (145)
Interestingly, these subtle sign changes preserve the basic properties of quaternions as rotation operators. Mathematically, the key consequence is the change of the sign of the cross-product in (13), which induces a change in the quaternion handedness (Shuster, 1993): Hamilton uses ij = k and is therefore right-handed, i.e., it turns vectors following the right-hand rule; JPL uses ji = k and is left-handed (Trawny and Roumeliotis, 2005). Being left-and right-handed rotations of opposite signs, we can say that their quaternions q lef t and q right are related by,

q left = q * right . (146)

Function of the rotation operator

We have seen how to rotate vectors in 3D. This is referred to in (Shuster, 1993) as the active interpretation, because operators (this affects all rotation operators) actively rotate vectors,

x = q active ⊗ x ⊗ q * active , x = R active x . (147)
Another way of seeing the effect of q and R over a vector x is to consider that the vector is steady but it is us who have rotated our point of view by an amount specified by q or R. This is called here frame transformation and it is referred to in (Shuster, 1993) as the passive interpretation, because vectors do not move,

x B = q passive ⊗ x A ⊗ q * passive , x B = R passive x A , (148)
where A and B are two Cartesian reference frames, and x A and x B are expressions of the same vector x in these frames. See further down for explanations and proper notations. The active and passive interpretations are governed by operators inverse of each other, that is,

q active = q * passive , R active = R passive .
Both Hamilton and JPL use the passive convention.

Direction cosine matrix A few authors understand the passive operator as not being a rotation operator, but rather an orientation specification, named the direction cosine matrix,

C =   c xx c xy c zx c xy c yy c zy c xz c yz c zz   , (149)
where each component c ij is the cosine of the angle between the axis i in the source frame and the axis j in the target frame. We have the identity,

C ≡ R passive . (150)

Direction of the rotation operator

In the passive case, a second source of interpretation is related to the direction in which the rotation matrix and quaternion operate, either converting from local to global frames, or from global to local.

Given two Cartesian frames G and L, we identify G and L as being the global and local frames. "Global" and "local" are relative definitions, i.e., G is global with respect to L, and L is local with respect to G -in other words, L is a frame specified in the reference frame G. 21 We specify q GL and R GL as being respectively the quaternion and rotation matrix transforming vectors from frame L to frame G, in the sense that a vector x L in frame L is expressed in frame G with the quaternion-and matrix-products

x G = q GL ⊗ x L ⊗ q * GL , x G = R GL x L . (151)
The opposite conversion, from G to L, is done with

x L = q LG ⊗ x G ⊗ q * LG , x L = R LG x G , (152)
where

q LG = q * GL , R LG = R GL . (153)
Hamilton uses local-to-global as the default specification of a frame L expressed in frame G,

q Hamilton q [with respect to][of ] = q [to][from] = q GL , (154)
while JPL uses the opposite, global-to-local conversion,

q JPL q [of ][with respect to] = q [to][from] = q LG . (155)
Notice that

q JPL q LG,lef t = q * LG,right = q GL,right q Hamilton , (156)
which is not particularly useful, but illustrates how easy it is to get confused when mixing conventions. Notice also that we can conclude that q JP L = q Hamilton , but this, far from being a beautiful result, is just the source of great confusion, because the equality is only present in the quaternion values, but the two quaternions, when employed in formulas, mean and represent different things.

4 Perturbations, derivatives and integrals 4.1 The additive and subtractive operators in SO(3)

In vector spaces R n , the addition and subtraction operations are performed with the regular sum '+' and minus '-' operations. In SO(3) this is not possible, but equivalent operators can be defined for establishing a proper calculus corpus.

We thus define the plus and minus operators, ⊕, , between elements R ∈ SO(3), and elements θ ∈ R 3 of the tangent space at R, as follows.

The plus operator. The 'plus' operator ⊕ : SO(3) × R 3 → SO(3) produces an element S of SO(3) which is the result of composing a reference element R of SO(3) with a (often small) rotation. This rotation is specified by a vector of θ ∈ R 3 in the vector space tangent to the SO(3) manifold at the reference element R, that is,

S = R ⊕ θ R • Exp(θ) R, S ∈ SO(3), θ ∈ R 3 . ( 157
)
Notice that this operator may be defined for any representation of SO(3). In particular, for the quaternion and rotation matrix we have,

q S = q R ⊕ θ = q R ⊗ Exp(θ) (158) R S = R R ⊕ θ = R R •Exp(θ) . (159)
The minus operator. The 'minus' operator : SO(3) × SO(3) → R 3 is the inverse of the above. It returns the vectorial angular difference θ ∈ R 3 between two elements of SO(3). This difference is expressed in the vector space tangent to the reference element R,

θ = S R Log(R -1 • S) R, S ∈ SO(3), θ ∈ R 3 , (160)
which for the quaternion and rotation matrix reads,

θ = q S q R = Log(q * R ⊗ q S ) (161) θ = R S R R = Log(R R R S ) . (162)
In both cases, notice that even though the vector difference θ is typically supposed to be small, the definitions above hold for any value of θ (up to the first coverage of the SO(3) manifold, that is, for angles θ < π).

4.2

The four possible derivative definitions 4.2.1 Functions from vector space to vector space

The scalar and vector cases follow the classical definition of the derivative: given a function f : R m → R n , we use {+, -} to define the derivative as

∂f (x) ∂x lim δx→0 f (x + δx) -f (x) δx ∈ R n×m (163)
Euler integration produces linear expressions of the form

f (x + ∆x) ≈ f (x) + ∂f (x) ∂x ∆x ∈ R n 4.2.2 Functions from SO(3) to SO(3)
Given a function f : SO(3) → SO(3) with R ∈ SO(3) and a local, small angular variation θ ∈ R 3 , we use {⊕, } to define the derivative as

∂f (R) ∂θ lim δθ→0 f (R ⊕ δθ) f (R) δθ ∈ R 3×3 (164) = lim δθ→0 Log f -1 (R) f (R Exp(δθ)) δθ (165)
Euler integration produces expressions of the form,

f (R ⊕ ∆θ) ≈ f (R) ⊕ ∂f (R) ∂θ ∆θ f (R) Exp ∂f (R) ∂θ ∆θ ∈ SO(3)
4.2.3 Functions from vector space to SO(3)

For the case of a function f : R m → SO(3), we use '+' for the vector perturbations, and ' ' for the SO(3) difference,

∂f (x) ∂x lim δx→0 f (x + δx) f (x) δx ∈ R 3×m (166) = lim δx→0 Log(f -1 (x)f (x + δx)) δx (167)
Euler integration produces expressions of the form,

f (x + ∆x) ≈ f (x) ⊕ ∂f (x) ∂x ∆x f (x) Exp ∂f (x) ∂x ∆x ∈ SO(3)

Functions from SO(3) to vector space

For the case of a function f : SO(3) → R n , we use '⊕' for the SO(3) perturbations, and '-' for the vector difference,

∂f (R) ∂θ lim δθ→0 f (R ⊕ δθ) -f (R) δθ ∈ R n×3 (168) = lim δθ→0 f (R Exp(δθ)) -f (R) δθ (169)
Euler integration produces expressions of the form,

f (R ⊕ δθ) ≈ f (R) + ∂f (R) ∂θ ∆θ f (R) + Exp ∂f (R) ∂θ ∆θ ∈ SO(3)

Useful, and very useful, Jacobians of the rotation

Let us consider a rotation to a vector a, of θ radians around the unit axis u. Let us express the rotation specification in three equivalent forms, namely θ = θu, q = q{θ} and R = R{θ}. We are interested in the Jacobians of the rotated result with respect to different magnitudes.

Jacobian with respect to the vector

The derivative of the rotation of a vector a with respect to this vector is trivial,

∂(q ⊗ a ⊗ q * ) ∂a = ∂(R a) ∂a = R . ( 170
)

Jacobian with respect to the quaternion

On the contrary, the derivative of the rotation with respect to the quaternion q is tricky.

For convenience, we use a lighter notation for the quaternion, q = [w v] = w + v. We make use of ( 35), ( 34), and the identity a

× (b × c) = (c × b) × a = (a c) b -(a b) c
, to develop the quaternion-based rotation (107) as follows,

a = q ⊗ a ⊗ q * = (w + v) ⊗ a ⊗ (w -v) = w 2 a + w(v ⊗ a -a ⊗ v) -v ⊗ a ⊗ v = w 2 a + 2w(v×a) -(-v a + v×a) ⊗ v = w 2 a + 2w(v×a) -(-v a) v + (v×a) ⊗ v = w 2 a + 2w(v×a) -(-v a) v - (v×a) v + (v×a)×v = w 2 a + 2w(v×a) -(-v a) v + (v v) a -(v a) v = w 2 a + 2w(v×a) + 2(v a) v -(v v) a . (171)
With this, we can extract the derivatives ∂a /∂w and ∂a /∂v, ∂a ∂w = 2(wa + v×a) (172)

∂a ∂v = -2w [a] × + 2(v a I + v a ) -2a v = 2(v a I + v a -a v -w [a] × ) , (173)
yielding

∂(q ⊗ a ⊗ q * ) ∂q = 2 w a + v×a v a I 3 + v a -a v -w [a] × ∈ R 3×4 . ( 174
) ✓ ✓ + ✓ R = exp(✓) exp(✓ + ✓) = exp(✓) exp( ) 1 = exp(0)

Right Jacobian of SO(3)

Let us consider (see Fig. 11) an element R ∈ SO(3) and a rotation vector θ ∈ R 3 such that R = Exp(θ). When θ is altered by an amount δθ, the element R varies. Expressing the variations of R in the tangent space of SO(3) at R with a rotation vector δφ ∈ R 3 , we have that (please see the figure, I am not inventing anything here)

Exp(θ) ⊕ δφ = Exp(θ + δθ) (175)
which might be written also as,

Exp(θ) • Exp(δφ) = Exp(θ + δθ) , (176)
and even

δφ = Log Exp(θ) -1 • Exp(θ + δθ) = Exp(θ + δθ) Exp(θ) . ( 177
)
In the limit, the variation of δφ as a function of δθ defines a Jacobian matrix

∂δφ ∂δθ = lim δθ→0 δφ δθ = lim δθ→0 Exp(θ + δθ) Exp(θ) δθ , ( 178
)
whose expression is a particular case of (166), that is, it is the derivative of the function f (θ) = Exp(θ), from R 3 to SO(3). This Jacobian matrix is known as the right Jacobian of SO(3), and is defined as,

J r (θ) ∂ Exp(θ) ∂θ . ( 179
)
Its expression is independent of the parametrization used, though it can indeed be expressed particularly for each parametrization. Using (166) we have,

J r (θ) = lim δθ→0 Exp(θ + δθ) Exp(θ) δθ (180) = lim δθ→0 Log(Exp(θ) Exp(θ + δθ)) δθ if using R (181) = lim δθ→0 Log(Exp(θ) * ⊗ Exp(θ + δθ)) δθ if using q . ( 182
)
The right Jacobian and its inverse can be computed in closed form (Chirikjian, 2012, page 40),

J r (θ) = I - 1 -cos θ θ 2 [θ] × + θ -sin θ θ 3 [θ] 2 × (183) J -1 r (θ) = I + 1 2 [θ] × + 1 θ 2 - 1 + cos θ 2 θ sin θ [θ] 2 × (184)
The right Jacobian of SO(3) has the following properties, for any θ and small δθ,

Exp(θ + δθ) ≈ Exp(θ) Exp(J r (θ)δθ) (185) Exp(θ) Exp(δθ) ≈ Exp(θ + J -1 r (θ) δθ) (186) Log(Exp(θ) Exp(δθ)) ≈ θ + J -1 r (θ) δθ (187)

Jacobian with respect to the rotation vector

The rotation of a vector a = R{θ} a with respect to the rotation vector θ is a function from R 3 to R 3 . Its derivative with respect to the rotation vector θ uses (163) and is developed from the previous result, using (185),

∂(q ⊗ a ⊗ q * ) ∂δθ = ∂(R a) ∂δθ = lim δθ→0 R{θ + δθ} a -R{θ} a δθ ← (163) = lim δθ→0 (R{θ} Exp(J r (θ) δθ) -R{θ})a δθ ← (185) = lim δθ→0 (R{θ}(I + [J r (θ) δθ] × ) -R{θ})a δθ = lim δθ→0 R{θ} [J r (θ) δθ] × a δθ = lim δθ→0 - R{θ} [a] × J r (θ) δθ δθ = -R{θ} [a] × J r (θ) ,
where R{θ} Exp(θ). Summarizing,

∂(q ⊗ a ⊗ q * ) ∂δθ = ∂(R a) ∂δθ = -R{θ} [a] × J r (θ) . ( 188
)
4.4 Perturbations, uncertainties, noise

Local perturbations

A perturbed orientation q may be expressed as the composition of the unperturbed orientation q with a small local perturbation ∆q L . Because of the Hamilton convention, this local perturbation appears at the right hand side of the composition product -we give also the matrix equivalent for comparison,

q = q ⊗ ∆q L , R = R ∆R L . ( 189
)
These local perturbation ∆q L (or ∆R L ) is easily obtained from its equivalent vector form ∆φ L = u∆φ L , defined in the tangent space, using the exponential map. This gives

qL = q L ⊗ Exp(∆φ L ) , RL = R L •Exp(∆φ L ) (190)
leading to an expression of the local perturbation

∆φ L = Log(q * L ⊗ qL ) = Log(R L • RL ) (191)
If the perturbation angle ∆φ L is small then the perturbation in quaternion and rotation matrix forms can be approximated by the Taylor expansions of ( 101) and ( 69) up to the linear terms,

∆q L ≈ 1 1 2 ∆φ L , ∆R L ≈ I + [∆φ L ] × . (192)
Perturbations can therefore be specified in the local vector space ∆φ L tangent to the SO(3) manifold at the actual orientation. It is convenient, for example, to express the covariances matrix of these perturbations in this vectorial space, that is, with a regular 3×3 covariance matrix.

Global perturbations

It is possible and indeed interesting to consider globally-defined perturbations, and likewise for the related derivatives. Global perturbations appear at the left hand side of the composition product, namely, 193) leading to an expression of the global perturbation

qG = Exp(∆φ G ) ⊗ q G , RG = Exp(∆φ G )•R G (
∆φ G = Log(q G ⊗ q * G ) = Log( RG •R G ) (194)
Again, these perturbations can be specified in the vector space ∆φ G tangent to the SO(3) manifold at the origin.

Time derivatives

Expressing the local perturbations in a vector space we can easily develop expressions for the time-derivatives. Just consider q = q(t) as the original state, q = q(t + ∆t) as the perturbed state, and apply the definition of the derivative

dq(t) dt lim ∆t→0 q(t + ∆t) -q(t) ∆t , (195)
to the above, with

ω L (t) dφ L (t) dt lim ∆t→0 ∆φ L ∆t , (196)
which, being ∆φ L a local angular perturbation, corresponds to the angular rates vector in the local frame defined by q.

The development of the time-derivative of the quaternion follows (an analogous reasoning would be used for the rotation matrix)

q lim ∆t→0 q(t + ∆t) -q(t) ∆t = lim ∆t→0 q ⊗ ∆q L -q ∆t = lim ∆t→0 q ⊗ 1 ∆φ L /2 - 1 0 ∆t = lim ∆t→0 q ⊗ 0 ∆φ L /2 ∆t = 1 2 q ⊗ 0 ω L . (197)
Defining

Ω(ω) [ω] R = 0 -ω ω -[ω] × =     0 -ω x -ω y -ω z ω x 0 ω z -ω y ω y -ω z 0 ω x ω z ω y -ω x 0     , (198)
we get from (197) and ( 17) (we give also its matrix equivalent)

q = 1 2 Ω(ω L ) q = 1 2 q ⊗ ω L , Ṙ = R [ω L ] × . (199)
These expressions are of course identical to (99) and ( 67), developed in the framework of the rotation group SO(3). Here, however, and interestingly, we are able to clearly refer the angular rate ω L to a particular reference frame, which in this case is the local frame defined by the orientation q or R. This has been possible now because we have given the operators q and R a precise geometrical meaning. From this viewpoint, (199) expresses the evolution of the orientation of a reference frame, when the angular rates are expressed locally in this frame.

The time-derivatives associated to global perturbations follow from a development analogous to (197), which results in

q = 1 2 ω G ⊗ q , Ṙ = [ω G ] × R , (200)
where

ω G (t) dφ G (t) dt (201)
is the angular rates vector expressed in the global frame. Eq. ( 200) expresses the evolution of the orientation of a reference frame, when the angular rates are expressed in the global reference frame.

Global-to-local relations

From the previous paragraph, it is worth noticing the following relation between local and global angular rates,

1 2 ω G ⊗ q = q = 1 2 q ⊗ ω L . (202)
Then, post-multiplying by the conjugate quaternion we have

ω G = q ⊗ ω L ⊗ q * = R ω L . (203)
Likewise, considering that ∆φ R ≈ ω∆t for small ∆t, we have that

∆φ G = q ⊗ ∆φ L ⊗ q * = R ∆φ L . (204)
That is, we can transform angular rates vectors ω and small angular perturbations ∆φ via frame transformation, using the quaternion or the rotation matrix, as if they were regular vectors. The same can be seen by posing ω = uω, or ∆φ = u∆φ, and noticing that the rotation axis vector u transforms normally, with

u G = q ⊗ u L ⊗ q * = R u L . (205)

Time-derivative of the quaternion product

We use the regular formula for the derivative of the product, (q

1 ⊗ q 2 ) = q1 ⊗ q 2 + q 1 ⊗ q2 , (R 1 R 2 ) = Ṙ1 R 2 + R 1 Ṙ2 , (206)
but noticing that, since the products are non commutative, we need to respect the order of the operands strictly. This means that (q 2 ) = 2 q ⊗ q , as it would be in the scalar case, but rather (q 2 ) = q ⊗ q + q ⊗ q .

(207)

Other useful expressions with the derivative

We can derive an expression for the local rotation rate

ω L = 2 q * ⊗ q , [ω L ] × = R Ṙ . ( 208
)
and the global rotation rate,

ω G = 2 q ⊗ q * , [ω G ] × = Ṙ R . (209)

Time-integration of rotation rates

Accumulating rotation over time in quaternion form is done by integrating the differential equation appropriate to the rotation rate definition, that is, (199) for a local rotation rate definition, and ( 200) for a global one. In the cases we are interested in, the angular rates are measured by local sensors, thus providing local measurements ω(t n ) at discrete times t n = n∆t. We concentrate here on this case only, for which we reproduce the differential equation (199),

q(t) = 1 2 q(t) ⊗ ω(t) . (210)
We develop zeroth-and first-order integration methods (Figs. 12 and 13), all based on the Taylor series of q(t n + ∆t) around the time t = t n . We note q q(t) and q n q(t n ), and the same for ω. The Taylor series reads,

q n+1 = q n + qn ∆t + 1 2! qn ∆t 2 + 1 3! ... q n ∆t 3 + 1 4! .... q n ∆t 4 + • • • . ( 211
)
The successive derivatives of q n above are easily obtained by repeatedly applying the expression of the quaternion derivative, (210), with ω = 0. We obtain

qn = 1 2 q n ω n (212a) qn = 1 2 2 q n ω 2 n + 1 2 q n ω (212b)
...

q n = 1 2 3 q n ω 3 n + 1 4 q n ω ω n + 1 2 q ω n ω (212c) q (i ≥ 4) n = 1 2 i q n ω i n + • • • , (212d)
Figure 12: Angular velocity approximations for the integral: Red: true velocity. Blue: zeroth order approximations (bottom to top: forward, midward and backward). Green: first order approximation.

! n ! n 1 q n 1 q n+1 q n ! n ! n+1 ! n 1 q n 1 q n+1 q n ! n ! n+1
q n 1 q n+1 q n where we have omitted the ⊗ signs for economy of notation, that is, all products and the powers of ω must be interpreted in terms of the quaternion product.

Zeroth order integration

Forward integration In the case where the angular rate ω n is held constant over the period [t n , t n+1 ], we have ω = 0 and (211) reduces to,

q n+1 = q n ⊗ 1 + 1 2 ω n ∆t + 1 2! 1 2 ω n ∆t 2 + 1 3! 1 2 ω n ∆t 3 + • • • , ( 213
)
where we identify the Taylor series (40) of the exponential e ωn∆t/2 . From (101), this exponential corresponds to the quaternion representing the incremental rotation ∆θ = ω n ∆t,

e ω∆t/2 = Exp(ω∆t) = q{ω∆t} = cos( ω ∆t/2) ω ω sin( ω ∆t/2)
, therefore,

q n+1 = q n ⊗ q{ω n ∆t} . ( 214
)
Backward integration We can also consider that the constant velocity over the period ∆t corresponds to ω n+1 , the velocity measured at the end of the period. This can be developed in a similar manner with a Taylor expansion of q n around t n+1 , leading to

q n+1 ≈ q n ⊗ q{ω n+1 ∆t} . ( 215
)
We want to remark here that this is the typical integration method when the arriving motion measurements are to be processed in real time, because the integration horizon corresponds to the last measurement (in this case, t n+1 , see Fig. 13). To make this more salient, we can re-label the time indices to use {n -1, n} instead of {n, n + 1}, and write,

q n = q n-1 ⊗ q{ω n ∆t} . ( 216
)
Midward integration Similarly, if the velocity is considered constant at the median rate over the period ∆t (which is not necessary the velocity at the midpoint of the period),

ω = ω n+1 + ω n 2 , (217)
we have,

q n+1 = q n ⊗ q{ ω∆t} . (218)

First order integration

The angular rate ω(t) is now linear with time. Its first derivative is constant, and all higher ones are zero,

ω = ω n+1 -ω n ∆t (219) ω = ... ω = • • • = 0 . ( 220
)
We can write the median rate ω in terms of ω n and ω,

ω = ω n + 1 2 ω∆t , (221)
and derive the expression of the powers of ω n appearing in the quaternion derivatives (212), in terms of the more convenient ω and ω,

ω n = ω - 1 2 ω∆t (222a) ω 2 n = ω 2 - 1 2 ω ω∆t - 1 2 ω ω∆t + 1 4 ω2 ∆t 2 (222b) ω 3 n = ω 3 - 3 2 ω 2 ω∆t + 3 4 ω ω2 ∆t 2 + 1 8 ω3 ∆t 3 (222c) ω 4 n = ω 4 + • • • . (222d)
Injecting them in the quaternion derivatives, and substituting in the Taylor series (211), we have after proper reordering,

q n+1 = q 1 + 1 2 ω∆t + 1 2! 1 2 ω∆t 2 + 1 3! 1 2 ω∆t 3 + • • • (223a) + q - 1 4 ω + 1 4 ω ∆t 2 (223b) + q - 1 16 ω ω - 1 16 ω ω + 1 24 ω ω + 1 12 ω ω ∆t 3 (223c) + q • • • ∆t 4 + • • • , ( 223d
)
where in (223a) we recognize the exponential series e ω∆t/2 = q{ ω∆t}, (223b) vanishes, and (223d) represents terms of high multiplicity that we are going to neglect. This yields after simplification (we recover now the normal ⊗ notation),

q n+1 = q n ⊗ q{ ω∆t} + ∆t 3 48 q n ⊗ ( ω ⊗ ω -ω ⊗ ω) + • • • . ( 224
)
Substituting ω and ω by their definitions ( 219) and ( 217) we get,

q n+1 = q n ⊗ q{ ω∆t} + ∆t 2 48 q n ⊗ ( ω n ⊗ ω n+1 -ω n+1 ⊗ ω n ) + • • • , ( 225
)
which is a result equivalent to (Trawny and Roumeliotis, 2005)'s, but using the Hamilton convention, and the quaternion product form instead of the matrix product form. Finally, since 34), we have the alternative form,

a v ⊗ b v -b v ⊗ a v = 2 a v ×b v , see (
q n+1 ≈ q n ⊗ q{ ω∆t} + ∆t 2 24 0 ω n × ω n+1 . ( 226
)
In this expression, the first term of the sum is the midward zeroth order integrator (218). The second term is a second-order correction that vanishes when ω n and ω n+1 are collinear,22 i.e., when the axis of rotation has not changed from t n to t n+1 .

Case of fixed rotation axis Let us write ω(t) = u(t) ω(t) and call u the axis of rotation.

In the case of a constant rotation axis u(t) = u, we have ω n × ω n+1 = 0 and therefore,

q n+1 = q n ⊗ q{u ω ∆t} . ( 227
)
This result is in fact interesting for cases not limited to first-order derivatives of ω(t). In effect, if the axis of rotation is constant, the infinitesimal contributions of rotation into the quaternion commute, i.e.,

exp(u ω 1 δt 1 ) exp(u ω 2 δt 2 ) = exp(u ω 2 δt 2 ) exp(u ω 1 δt 1 ) = exp(u (ω 1 δt 1 + ω 2 δt 2 )) ,
and thus we have the identity,

q n+1 = q n ⊗ exp u 2 t n+1 tn ω(t) δt (228a) = q n ⊗ exp(u ∆θ n /2) (228b) = q n ⊗ q{u ∆θ n } . ( 228c
)
with ∆θ n = t n+1 tn ω(t)dt ∈ R the total angle rotated during the interval [t n , t n+1 ].

Case of varying rotation axis Clearly, the second term of the sum in ( 226) captures through ω n × ω n+1 = 0 the effect that a varying rotation axis has on the integrated orientation. For its practical usage, we notice that given usual IMU sampling times ∆t ≤ 0.01s, and the usual near-collinearity of ω n and ω n+1 due to inertia, this second-order term takes values of the order of 10 -6 ω 2 , or easily smaller. Terms with higher multiplicities of ω∆t are even smaller and have been neglected. Please note also that, while all zeroth-order integrators result in unit quaternions by construction (because they are computed as the product of two unit quaternions), this is not the case for the first-order integrator due to the sum in (226). Hence, when using the first-order integrator, and even if the summed term is small as stated, users should take care to check the evolution of the quaternion norm over time, and eventually re-normalize the quaternion if needed, using quaternion updates of the form q ← q/ q . Only if the constant axis assumption holds, then (227) holds too and this normalization is no longer necessary.

5 Error-state kinematics for IMU-driven systems

Motivation

We wish to write the error-estate equations of the kinematics of an inertial system integrating accelerometer and gyrometer readings with bias and noise, using the Hamilton quaternion to represent the orientation in space or attitude.

Accelerometer and gyrometer readings come typically from an Inertial Measurement Unit (IMU). Integrating IMU readings leads to dead-reckoning positioning systems, which drift with time. Avoiding drift is a matter of fusing this information with absolute position readings such as GPS or vision.

The error-state Kalman filter (ESKF) is one of the tools we may use for this purpose. Within the Kalman filtering paradigm, these are the most remarkable assets of the ESKF (Madyastha et al., 2011):

• The orientation error-state is minimal (i.e., it has the same number of parameters as degrees of freedom), avoiding issues related to over-parametrization (or redundancy) and the consequent risk of singularity of the involved covariances matrices, resulting typically from enforcing constraints.

• The error-state system is always operating close to the origin, and therefore far from possible parameter singularities, gimbal lock issues, or the like, providing a guarantee that the linearization validity holds at all times.

• The error-state is always small, meaning that all second-order products are negligible. This makes the computation of Jacobians very easy and fast. Some Jacobians may even be constant or equal to available state magnitudes.

• The error dynamics are slow because all the large-signal dynamics have been integrated in the nominal-state. This means that we can apply KF corrections (which are the only means to observe the errors) at a lower rate than the predictions.

The error-state Kalman filter explained

In error-state filter formulations, we speak of true-, nominal-and error-state values, the true-state being expressed as a suitable composition (linear sum, quaternion product or matrix product) of the nominal-and the error-states. The idea is to consider the nominalstate as large-signal (integrable in non-linear fashion) and the error-state as small signal (thus linearly integrable and suitable for linear-Gaussian filtering).

The error-state filter can be explained as follows. On one side, high-frequency IMU data u m is integrated into a nominal-state x. This nominal state does not take into account the noise terms w and other possible model imperfections. As a consequence, it will accumulate errors. These errors are collected in the error-state δx and estimated with the Error-State Kalman Filter (ESKF), this time incorporating all the noise and perturbations. The errorstate consists of small-signal magnitudes, and its evolution function is correctly defined by a (time-variant) linear dynamic system, with its dynamic, control and measurement matrices computed from the values of the nominal-state. In parallel with integration of the nominalstate, the ESKF predicts a Gaussian estimate of the error-state. It only predicts, because by now no other measurement is available to correct these estimates. The filter correction is performed at the arrival of information other than IMU (e.g. GPS, vision, etc.), which is able to render the errors observable and which happens generally at a much lower rate than the integration phase. This correction provides a posterior Gaussian estimate of the error-state. After this, the error-state's mean is injected into the nominal-state, then reset to zero. The error-state's covariances matrix is conveniently updated to reflect this reset. The system goes on like this forever.

Velocity v t v δv v t = v + δv Quaternion ( 2,3 ) q t q δq q t = q ⊗ δq Rotation matrix ( 2,3 ) R t R δR R t = R δR Angles vector ( 4 ) δθ δq = e δθ/2 δR = e [δθ] × Accelerometer bias a bt a b δa b a bt = a b + δa b a w Gyrometer bias ω bt ω b δω b ω bt = ω b + δω b ω w Gravity vector g t g δg g t = g + δg Acceleration a t a m a n Angular rate ω t ω m ω n
( 1 ) the symbol ⊕ indicates a generic composition

( 2 ) indicates non-minimal representations

( 3 ) see Table 4 for the composition formula in case of globally-defined angular errors ( 4 ) exponentials defined as in ( 101) and (69, 77)

System kinematics in continuous time

The definition of all the involved variables is summarized in Table 3. Two important decisions regarding conventions are worth mentioning:

• The angular rates ω are defined locally with respect to the nominal quaternion. This allows us to use the gyrometer measurements ω m directly, as they provide bodyreferenced angular rates.

• The angular error δθ is also defined locally with respect to the nominal orientation. This is not necessarily the optimal way to proceed, but it corresponds to the choice in most IMU-integration works -what we could call the classical approach. There exists evidence (Li and Mourikis, 2012) that a globally-defined angular error has better properties. This will be explored too in the present document, Section 7, but most of the developments, examples and algorithms here are based in this locallydefined angular error.

which we may name ẋt = f t (x t , u, w). This system has state x t , is governed by IMU noisy readings u m , and is perturbed by white Gaussian noise w, all defined by

x t =         p t v t q t a bt ω bt g t         u = a m -a n ω m -ω n w = a w ω w . ( 235
)
It is to note in the above formulation that the gravity vector g t is going to be estimated by the filter. It has a constant evolution equation, (234f), as corresponds to a magnitude that is known to be constant. The system starts at a fixed and arbitrarily known initial orientation q t (t = 0) = q 0 , which, being generally not in the horizontal plane, makes the initial gravity vector generally unknown. For simplicity it is usually taken q 0 = (1, 0, 0, 0) and thus R 0 = R{q 0 } = I. We estimate g t expressed in frame q 0 , and not q t expressed in a horizontal frame, so that the initial uncertainty in orientation is transferred to an initial uncertainty on the gravity direction. We do so to improve linearity: indeed, equation ( 234b) is now linear in g, which carries all the uncertainty, and the initial orientation q 0 is known without uncertainty, so that q starts with no uncertainty. Once the gravity vector is estimated the horizontal plane can be recovered and, if desired, the whole state and recovered motion trajectories can be re-oriented to reflect the estimated horizontal. See (Lupton and Sukkarieh, 2009) for further justification. This is of course optional, and the reader is free to remove all equations related to graviy from the system and adopt a more classical approach of considering g (0, 0, -9.8xx), with xx the appropriate decimal digits of the gravity vector on the site of the experiment, and an uncertain initial orientation q 0 .

The nominal-state kinematics

The nominal-state kinematics corresponds to the modeled system without noises or perturbations,

ṗ = v (236a) v = R(a m -a b ) + g (236b) q = 1 2 q ⊗ (ω m -ω b ) (236c) ȧb = 0 (236d) ωb = 0 (236e) ġ = 0. (236f)

The error-state kinematics

The goal is to determine the linearized dynamics of the error-state. For each state equation, we write its composition (in Table 3), solving for the error state and simplifying all second-order infinitesimals. We give here the full error-state dynamic system and proceed afterwards with comments and proofs.

δ p = δv (237a) δ v = -R [a m -a b ] × δθ -Rδa b + δg -Ra n (237b) δ θ = -[ω m -ω b ] × δθ -δω b -ω n (237c) δ a b = a w (237d) δω b = ω w (237e) δ g = 0. (237f)
Equations ( 237a), ( 237d), ( 237e) and (237f), respectively of position, both biases, and gravity errors, are derived from linear equations and their error-state dynamics is trivial.

As an example, consider the true and nominal position equations ( 234a) and (236a), their composition p t = p + δp from Table 3, and solve for δ p to obtain (237a). Equations ( 237b) and (237c), of velocity and orientation errors, require some nontrivial manipulations of the non-linear equations ( 234b) and (234c) to obtain the linearized dynamics. Their proofs are developed in the following two sections.

Equation (237b):

The linear velocity error. We wish to determine δ v, the dynamics of the velocity errors. We start with the following relations

R t = R(I + [δθ] × ) + O( δθ 2 ) (238) v = Ra B + g, (239)
where ( 238) is the small-signal approximation of R t , and in (239) we rewrote (236b) but introducing a B and δa B , defined as the large-and small-signal accelerations in body frame,

a B a m -a b (240) δa B -δa b -a n (241)
so that we can write the true acceleration in inertial frame as a composition of large-and small-signal terms,

a t = R t (a B + δa B ) + g + δg. (242)
We proceed by writing the expression (234b) of vt in two different forms (left and right developments), where the terms O( δθ 2 ) have been ignored,

v + δ v = vt = R(I + [δθ] × )(a B + δa B ) + g + δg Ra B + g + δ v = = Ra B + Rδa B + R [δθ] × a B + R [δθ] × δa B + g + δg
This leads after removing Ra B + g from left and right to

δ v = R(δa B + [δθ] × a B ) + R [δθ] × δa B + δg (243)
so that ω t can be written with a nominal part and an error part,

ω t = ω + δω. (254)
We proceed by computing qt by two different means (left and right developments)

(q ⊗ δq) = qt = 1 2 q t ⊗ ω t q ⊗ δq + q ⊗ δ q = = 1 2 q ⊗ δq ⊗ ω t 1 2 q ⊗ ω ⊗ δq + q ⊗ δ q =
simplifying the leading q and isolating δ q we obtain

0 δ θ = 2 δ q = δq ⊗ ω t -ω ⊗ δq = [q] R (ω t )δq -[q] L (ω)δq = 0 -(ω t -ω) (ω t -ω) -[ω t + ω] × 1 δθ/2 + O( δθ 2 ) = 0 -δω δω -[2ω + δω] × 1 δθ/2 + O( δθ 2 ) (255)
which results in one scalar-and one vector-equalities

0 = δω δθ + O(|δθ| 2 ) (256a) δ θ = δω -[ω] × δθ - 1 2 [δω] × δθ + O( δθ 2 ). (256b)
The first equation leads to δω δθ = O( δθ 2 ), which is formed by second-order infinitesimals, not very useful. The second equation yields, after neglecting all second-order terms,

δ θ = -[ω] × δθ + δω (257)
and finally, recalling ( 252) and ( 253), we get the linearized dynamics of the angular error,

δ θ = -[ω m -ω b ] × δθ -δω b -ω n . (258)

System kinematics in discrete time

The differential equations above need to be integrated into differences equations to account for discrete time intervals ∆t > 0. The integration methods may vary. In some cases, one will be able to use exact closed-form solutions. In other cases, numerical integration of varying degree of accuracy may be employed. Please refer to the Appendices for pertinent details on integration methods. Integration needs to be done for the following sub-systems:

1. The nominal state.

2. The error-state.

(a) The deterministic part: state dynamics and control.

(b) The stochastic part: noise and perturbations.

The nominal state kinematics

We can write the differences equations of the nominal-state as

p ← p + v ∆t + 1 2 (R(a m -a b ) + g) ∆t 2 (259a) v ← v + (R(a m -a b ) + g) ∆t (259b) q ← q ⊗ q{(ω m -ω b ) ∆t} (259c) a b ← a b (259d) ω b ← ω b (259e) g ← g , (259f)
where x ← f (x, •) stands for a time update of the type

x k+1 = f (x k , • k ), R
R{q} is the rotation matrix associated to the current nominal orientation q, and q{v} is the quaternion associated to the rotation v, according to (101).

We can also use more precise integration, please see the Appendices for more information.

The error-state kinematics

The deterministic part is integrated normally (in this case we follow the methods in App. C.2), and the integration of the stochastic part results in random impulses (see App. E), thus,

δp ← δp + δv ∆t (260a) δv ← δv + (-R [a m -a b ] × δθ -Rδa b + δg)∆t + v i (260b) δθ ← R {(ω m -ω b )∆t}δθ -δω b ∆t + θ i (260c) δa b ← δa b + a i (260d) δω b ← δω b + ω i (260e) δg ← δg . (260f)
Here, v i , θ i , a i and ω i are the random impulses applied to the velocity, orientation and bias estimates, modeled by white Gaussian processes. Their mean is zero, and their covariances matrices are obtained by integrating the covariances of a n , ω n , a w and ω w over the step time ∆t (see App. E),

V i = σ 2 ãn ∆t 2 I [m 2 /s 2 ] (261) Θ i = σ 2 ωn ∆t 2 I [rad 2 ] (262) A i = σ 2 aw ∆tI [m 2 /s 4 ] (263) Ω i = σ 2 ωw ∆tI [rad 2 /s 2 ] ( 264
)
where

σ ãn [m/s 2 ], σ ωn [rad/s], σ aw [m/s 2 √ s] and σ ωw [rad/s √ s]
are to be determined from the information in the IMU datasheet, or from experimental measurements.

The error-state Jacobian and perturbation matrices

The Jacobians are obtained by simple inspection of the error-state differences equations in the previous section.

To write these equations in compact form, we consider the nominal state vector x, the error state vector δx, the input vector u m , and the perturbation impulses vector i, as follows (see App. E.1 for details and justifications),

x =         p v q a b ω b g         , δx =         δp δv δθ δa b δω b δg         , u m = a m ω m , i =     v i θ i a i ω i     (265)
The error-state system is now

δx ← f (x, δx, u m , i) = F x (x, u m )•δx + F i •i, (266)
The ESKF prediction equations are written: δx ← F x (x, u m )• δx (267)

P ← F x P F x + F i Q i F i , (268)
where δx ∼ N { δx, P}25 ; F x and F i are the Jacobians of f () with respect to the error and perturbation vectors; and Q i is the covariances matrix of the perturbation impulses.

The expressions of the Jacobian and covariances matrices above are detailed below. All state-related values appearing herein are extracted directly from the nominal state.

F x = ∂f ∂δx x,um =         I I∆t 0 0 0 0 0 I -R [a m -a b ] × ∆t -R∆t 0 I∆t 0 0 R {(ω m -ω b )∆t} 0 -I∆t 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I         (269)
F i = ∂f ∂i x,um =         0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0         , Q i =     V i 0 0 0 0 Θ i 0 0 0 0 A i 0 0 0 0 Ω i     . (270)
Please note particularly that F x is the system's transition matrix, which can be approximated to different levels of precision in a number of ways. We showed here one of its simplest forms (the Euler form). Se Appendices B to D for further reference.

Please note also that, being the mean of the error δx initialized to zero, the linear equation ( 267) always returns zero. You should of course skip line (267) in your code. I recommend that you write it, though, but that you comment it out so that you are sure you did not forget anything.

And please note, finally, that you should NOT skip the covariance prediction (268)!! In effect, the term F i Q i F i is not null and therefore this covariance grows continuously -as it must be in any prediction step.

Fusing IMU with complementary sensory data

At the arrival of other kind of information than IMU, such as GPS or vision, we proceed to correct the ESKF. In a well-designed system, this should render the IMU biases observable and allow the ESKF to correctly estimate them. There are a myriad of possibilities, the most popular ones being GPS + IMU, monocular vision + IMU, and stereo vision + IMU. In recent years, the combination of visual sensors with IMU has attracted a lot of attention, and thus generated a lot of scientific activity. These vision + IMU setups are very interesting for use in GPS-denied environments, and can be implemented on mobile devices (typically smart phones), but also on UAVs and other small, agile platforms.

While the IMU information has served so far to make predictions to the ESKF, this other information is used to correct the filter, and thus observe the IMU bias errors. The correction consists of three steps:

1. observation of the error-state via filter correction, 2. injection of the observed errors into the nominal state, and 3. reset of the error-state.

These steps are developed in the following sections.

Observation of the error state via filter correction

Suppose as usual that we have a sensor that delivers information that depends on the state, such as

y = h(x t ) + v , (271)
where h() is a general nonlinear function of the system state (the true state), and v is a white Gaussian noise with covariance V,

v ∼ N {0, V} . ( 272
)
Our filter is estimating the error state, and therefore the filter correction equations26 ,

K = PH (HPH + V) -1 (273) δx ← K(y -h(x t )) (274) P ← (I -KH)P (275)
require the Jacobian matrix H to be defined with respect to the error state δx, and evaluated at the best true-state estimate xt = x ⊕ δx. As the error state mean is zero at this stage (we have not observed it yet), we have xt = x and we can use the nominal error x as the evaluation point, leading to

H ≡ ∂h ∂δx x . (276)

Jacobian computation for the filter correction

The Jacobian above might be computed in a number of ways. The most illustrative one is by making use of the chain rule,

H ∂h ∂δx x = ∂h ∂x t x ∂x t ∂δx x = H x X δx . (277)
Here, H x ∂h ∂xt x is the standard Jacobian of h() with respect to its own argument (i.e., the Jacobian one would use in a regular EKF). This first part of the chain rule depends on the measurement function of the particular sensor used, and is not presented here.

The second part, X δx ∂xt ∂δx x , is the Jacobian of the true state with respect to the error state. This part can be derived here as it only depends on the ESKF composition of states. We have the derivatives,

X δx =           ∂(p+δp) ∂δp ∂(v+δv) ∂δv 0 ∂(q⊗δq) ∂δθ ∂(a b +δa b ) ∂δa b 0 ∂(ω b +δω b ) ∂δω b ∂(g+δg) ∂δg           (278)
which results in all identity 3 × 3 blocks (for example, ∂(p+δp) ∂δp = I 3 ) except for the 4 × 3 quaternion term Q δθ = ∂(q ⊗ δq)/∂δθ. Therefore we have the form,

X δx ∂x t ∂δx x =   I 6 0 0 0 Q δθ 0 0 0 I 9   (279)
Using the chain rule, equations (17-19), and the limit δq -→

δθ→0 1 1 2 δθ
, the quaternion term Q δθ may be derived as follows,

Q δθ ∂(q ⊗ δq) ∂δθ q = ∂(q ⊗ δq) ∂δq q ∂δq ∂δθ δθ=0 = ∂([q] L δq) ∂δq q ∂ 1 1 2 δθ ∂δθ δθ=0 = [q] L 1 2     0 0 0 1 0 0 0 1 0 0 0 1     ,
which leads to

Q δθ = 1 2    
-q x -q y -q z q w -q z q y q z q w -q x -q y q x q w     .

(280)

Injection of the observed error into the nominal state

After the ESKF update, the nominal state gets updated with the observed error state using the appropriate compositions (sums or quaternion products, see Table 3),

x ← x ⊕ δx , (281)
that is,

p ← p + δp (282a) v ← v + δv (282b) q ← q ⊗ q{ δθ} (282c) a b ← a b + δ a b (282d) ω b ← ω b + δ ω b (282e)
g ← g + δg (282f)

ESKF reset

After error injection into the nominal state, the error state mean δx gets reset. This is especially relevant for the orientation part, as the new orientation error will be expressed locally with respect to the orientation frame of the new nominal state. To make the ESKF update complete, the covariance of the error needs to be updated according to this modification.

Let us call the error reset function g(). It is written as follows,

δx ← g(δx) = δx δx , (283)
where stands for the composition inverse of ⊕. The ESKF error reset operation is thus, δx ← 0 (284)

P ← G P G . ( 285
)
where G is the Jacobian matrix defined by,

G ∂g ∂δx δx . (286)
Similarly to what happened with the update Jacobian above, this Jacobian is the identity on all diagonal blocks except in the orientation error. We give here the full expression and proceed in the following section with the derivation of the orientation error block,

∂δθ + /∂δθ = I -1 2 δθ × , G =    I 6 0 0 0 I -1 2 δθ × 0 0 0 I 9    . (287)
In major cases, the error term δθ can be neglected, leading simply to a Jacobian G = I 18 , and thus to a trivial error reset. This is what most implementations of the ESKF do. The expression here provided should produce more precise results, which might be of interest for reducing long-term error drift in odometry systems.

Jacobian of the reset operation with respect to the orientation error

We want to obtain the expression of the new angular error δθ + with respect to the old error δθ and the observed error δθ. Consider these facts:

• The true orientation does not change on error reset, i.e., q + t = q t . This gives:

q + ⊗ δq + = q ⊗ δq . (288)
• The observed error mean has been injected into the nominal state (see ( 282c) and ( 123)):

q + = q ⊗ δq . ( 289
)
Combining both identities we obtain an expression of δq + ,

δq + = (q + ) * ⊗ q ⊗ δq = (q ⊗ δq) * ⊗ q ⊗ δq = δq * ⊗ δq = [ δq * ] L • δq . (290)
Considering that δq * ≈ 1 -1 2 δθ , the identity above can be expanded as

1 1 2 δθ + =   1 1 2 δθ -1 2 δθ I -1 2 δθ ×   • 1 1 2 δθ + O( δθ 2 ) , (291)
which gives one scalar-and one vector-equations,

1 4 δθ δθ = O( δθ 2 ) (292a) δθ + = -δθ + I - 1 2 δθ × δθ + O( δθ 2 ) , (292b)
among which the first one is not very informative in that it is only a relation of infinitesimals. One can show from the second equation that δθ + = 0, which is what we expect from the reset operation. The Jacobian is obtained by simple inspection,

∂δθ + ∂δθ = I - 1 2 δθ × . ( 293
)
7 The ESKF using global angular errors

We explore in this section the implications of having the angular error defined in the global reference, as opposed to the local definition we have used so far. We retrace the development of sections 5 and 6, and particularize the subsections that present changes with respect to the new definition.

A global definition of the angular error δθ implies a composition on the left hand side, i.e., q t = δq ⊗ q = q{δθ} ⊗ q .

We remark for the sake of completeness that we keep the local definition of the angular rates vector ω, i.e., q = 1 2 q ⊗ ω in continuous time, and therefore q ← q ⊗ q{ω∆t} in discrete time, regardless of the angular error being defined globally. This is so for convenience, as the measure of the angular rates provided by the gyrometers is in body frame, that is, local. and finally, recalling ( 240) and ( 241) and rearranging, we obtain the expression of the derivative of the velocity error,

δ v = -[R(a m -a b )] × δθ -Rδa b + δg -Ra n (299)
Equation ( 294c): The orientation error. We start by writing the true-and nominaldefinitions of the quaternion derivatives,

qt = 1 2 q t ⊗ ω t (300) q = 1 2 q ⊗ ω , (301)
and remind that we are using a globally-defined angular error, i.e.,

q t = δq ⊗ q . ( 302
)
As we did for the locally-defined error case, we also group large-and small-signal angular rates (252-253). We proceed by computing qt by two different means (left and right developments),

(δq ⊗ q) = qt = 1 2 q t ⊗ ω t δ q ⊗ q + δq ⊗ q = = 1 2 δq ⊗ q ⊗ ω t δ q ⊗ q + 1 2 δq ⊗ q ⊗ ω =
Having ω t = ω + δω, this reduces to

δ q ⊗ q = 1 2 δq ⊗ q ⊗ δω . (303)
Right-multiplying left and right terms by q * , and recalling that q ⊗ δω ⊗ q * ≡ Rδω, we can further develop as follows,

δ q = 1 2 δq ⊗ q ⊗ δω ⊗ q * = 1 2 δq ⊗ (Rδω) = 1 2 δq ⊗ δω G , (304)
with δω G Rδω the small-signal angular rate expressed in the global frame. Then,

0 δ θ = 2 δ q = δq ⊗ δω G = Ω(δω G ) δq = 0 -δω G δω G -[δω G ] × 1 δθ/2 + O( δθ 2 ) , (305)
which results in one scalar-and one vector-equalities

0 = δω G δθ + O(|δθ| 2 ) (306a) δ θ = δω G - 1 2 [δω G ] × δθ + O( δθ 2 ). ( 306b
)
The first equation leads to δω G δθ = O( δθ 2 ), which is formed by second-order infinitesimals, not very useful. The second equation yields, after neglecting all second-order terms,

δ θ = δω G = Rδω . (307)
Finally, recalling (253), we obtain the linearized dynamics of the global angular error,

δ θ = -Rδω b -Rω n . ( 308
)
7.2 System kinematics in discrete time

The nominal state

The nominal state equations do not involve errors and are therefore the same as in the case where the orientation error is defined locally.

The error state

Using Euler integration, we obtain the following set of differences equations,

δp ← δp + δv ∆t (309a) δv ← δv + (-[R(a m -a b )] × δθ -Rδa b + δg)∆t + v i (309b) δθ ← δθ -Rδω b ∆t + θ i (309c) δa b ← δa b + a i (309d) δω b ← δω b + ω i (309e) δg ← δg.
(309f)

The error state Jacobian and perturbation matrices

The Transition matrix is obtained by simple inspection of the equations above,

F x =            I I∆t 0 0 0 0 0 I -[R(a m -a b )] × ∆t -R∆t 0 I∆t 0 0 I 0 -R∆t 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I            . ( 310
)
We observe three changes with respect to the case with a locally-defined angular error (compare the boxed terms in the Jacobian above to the ones in ( 269)); these changes are summarized in Table 4.

The perturbation Jacobian and the perturbation matrix are unchanged after considering isotropic noises and the developments of App. E,

F i =         0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0         , Q i =     V i 0 0 0 0 Θ i 0 0 0 0 A i 0 0 0 0 Ω i     .
(311)

Fusing with complementary sensory data

The fusing equations involving the ESKF machinery vary only slightly when considering global angular errors. We revise these variations in the error state observation via ESKF correction, the injection of the error into the nominal state, and the reset step.

Error state observation

The only difference with respect to the local error definition is in the Jacobian block of the observation function that relates the orientation to the angular error. This new block is developed below.

Using (17-19) and the first-order approximation δq → 1 1 2 δθ , the quaternion term Q δθ may be derived as follows,

Q δθ ∂(δq ⊗ q) ∂δθ q = ∂(δq ⊗ q) ∂δq q ∂δq ∂δθ δθ=0 (312a) = [q] R 1 2     0 0 0 1 0 0 0 1 0 0 0 1     (312b) = 1 2    
-q x -q y -q z q w q z -q y -q z q w q x q y -q x q w     .

(312c)

Injection of the observed error into the nominal state

The composition x ← x ⊕ δx of the nominal and error states is depicted as follows,

p ← p + δp (313a) v ← v + δv (313b) q ← q{ δθ} ⊗ q (313c) a b ← a b + δa b (313d) ω b ← ω b + δω b (313e) g ← g + δg . ( 313f
)
where only the equation for the quaternion update has been affected. This is summarized in Table 4.

ESKF reset

The ESKF error mean is reset, and the covariance updated, according to, δx ← 0 (314)

P ← GPG (315)
with the Jacobian

G =    I 6 0 0 0 I + 1 2 δθ × 0 0 0 I 9    (316)
whose non-trivial term is developed as follows. Our goal is to obtain the expression of the new angular error δθ + with respect to the old error δθ. We consider these facts:

• The true orientation does not change on error reset, i.e., q + t ≡ q t . This gives: δq + ⊗ q + = δq ⊗ q .

(317)

• The observed error mean has been injected into the nominal state (see ( 282c) and ( 123)):

q + = δq ⊗ q . ( 318
)
Combining both identities we obtain an expression of the new orientation error with respect to the old one and the observed error δq,

δq + = δq ⊗ δq * = [ δq * ] R • δq . ( 319
)
Considering that δq * ≈ 1 -1 2 δθ , the identity above can be expanded as

1 1 2 δθ + =   1 1 2 δθ -1 2 δθ I + 1 2 δθ ×   • 1 1 2 δθ + O( δθ 2 ) (320)
t q t = q ⊗ δq q t = δq ⊗ q Euler integration ∂δv + /∂δθ -R [a m -a b ] × ∆t -[R(a m -a b )] × ∆t ∂δθ + /∂δθ R {(ω m -ω b )∆t} I ∂δθ + /∂δω b -I∆t -R∆t Error observation Q δθ 1 2          -q x -q y -q z q w -q z q y q z q w -q x -q y q x q w          1 2         
-q x -q y -q z q w q z -q y -q z q w q x q y -q

x q w          Error injection q ← q ⊗ q{ δθ} q ← q{ δθ} ⊗ q Error reset ∂δθ + /∂δθ I -1 2 δθ × I + 1 2 δθ
× which gives one scalar-and one vector-equations,

1 4 δθ δθ = O( δθ 2 ) (321a) δθ + = -δθ + I + 1 2 δθ × δθ + O( δθ 2 ) (321b)
among which the first one is not very informative in that it is only a relation of infinitesimals. One can show from the second equation that δθ + = 0, which is what we expect from the reset operation. The Jacobian is obtained by simple inspection,

∂δθ + ∂δθ = I + 1 2 δθ × . (322)
The difference with respect to the local error case is summarized in Table 4.

A Runge-Kutta numerical integration methods

We aim at integrating nonlinear differential equations of the form

ẋ = f (t, x) (323)
over a limited time interval ∆t, in order to convert them to a differences equation, i.e.,

x(t + ∆t) = x(t) + t+∆t t f (τ, x(τ ))dτ , (324)
or equivalently, if we assume that t n = n∆t and x n x(t n ),

x n+1 = x n + (n+1)∆t n∆t f (τ, x(τ ))dτ . (325)
One of the most utilized family of methods is the Runge-Kutta methods (from now on, RK). These methods use several iterations to estimate the derivative over the interval, and then use this derivative to integrate over the step ∆t.

In the sections that follow, several RK methods are presented, from the simplest one to the most general one, and are named according to their most common name. NOTE: All the material here is taken from the Runge-Kutta method entry in the English Wikipedia.

A.1 The Euler method

The Euler method assumes that the derivative f (•) is constant over the interval, and therefore

x n+1 = x n + ∆t•f (t n , x n ) . (326)
Put as a general RK method, this corresponds to a single-stage method, which can be depicted as follows. Compute the derivative at the initial point,

k 1 = f (t n , x n ) , (327)
and use it to compute the integrated value at the end point,

x n+1 = x n + ∆t•k 1 . (328)

A.2 The midpoint method

The midpoint method assumes that the derivative is the one at the midpoint of the interval, and makes one iteration to compute the value of x at this midpoint, i.e.,

x n+1 = x n + ∆t•f t n + 1 2 ∆t , x n + 1 2 ∆t•f (t n , x n ) . (329)
The midpoint method can be explained as a two-step method as follows. First, use the Euler method to integrate until the midpoint, using k 1 as defined previously,

k 1 = f (t n , x n ) (330) x(t n + 1 2 ∆t) = x n + 1 2 ∆t•k 1 . (331)
Then use this value to evaluate the derivative at the midpoint, k 2 , leading to the integration

k 2 = f (t n + 1 2 ∆t , x(t n + 1 2 ∆t)) (332) x n+1 = x n + ∆t•k 2 .
(333)

A.3 The RK4 method

This is usually referred to as simply the Runge-Kutta method. It assumes evaluation values for f () at the start, midpoint and end of the interval. And it uses four stages or iterations to compute the integral, with four derivatives, k 1 . . . k 4 , that are obtained sequentially. These derivatives, or slopes, are then weight-averaged to obtain the 4th-order estimate of the derivative in the interval. The RK4 method is better specified as a small algorithm than a one-step formula like the two methods above. The RK4 integration step is,

x n+1 = x n + ∆t 6 k 1 + 2k 2 + 2k 3 + k 4 , (334)
that is, the increment is computed by assuming a slope which is the weighted average of the slopes k 1 , k 2 , k 3 , k 4 , with

k 1 = f (t n , x n ) (335) k 2 = f t n + 1 2 ∆t , x n + ∆t 2 k 1 (336) k 3 = f t n + 1 2 ∆t , x n + ∆t 2 k 2 (337) k 4 = f t n + ∆t , x n + ∆t•k 3 . (338)
The different slopes have the following interpretation:

• k 1 is the slope at the beginning of the interval, using x n , (Euler's method);

• k 2 is the slope at the midpoint of the interval, using x n + 1 2 ∆t•k 1 , (midpoint method); • k 3 is again the slope at the midpoint, but now using x n + 1 2 ∆t•k 2 ; • k 4 is the slope at the end of the interval, using x n + ∆t•k 3 .

A.4 General Runge-Kutta method

More elaborated RK methods are possible. They aim at either reduce the error and/or increase stability. They take the general form

x n+1 = x n + ∆t s i=1 b i k i , (339)

B.1 Integration of the angular error

For example, consider the angular error dynamics without bias and noise (a cleaned version of Eq. ( 237c)),

δ θ = -[ω] × δθ (344)
Its transition matrix can be written as a Taylor series,

Φ = e -[ω] × ∆t (345) = I -[ω] × ∆t + 1 2 [ω] 2 × ∆t 2 - 1 3! [ω] 3 × ∆t 3 + 1 4! [ω] 4 × ∆t 4 -. . . (346)
Now defining ω∆t u∆θ, the unitary axis of rotation and the rotated angle, and applying (75), we can group terms and get

Φ = I -[u] × ∆θ + 1 2 [u] 2 × ∆θ 2 - 1 3! [u] 3 × ∆θ 3 + 1 4! [y] 4 × ∆θ 4 -. . . = I -[u] × ∆θ - ∆θ 3 3! + ∆θ 5 5! -• • • + [u] 2 × ∆θ 2 2! - ∆θ 4 4! + ∆θ 6 6! -• • • = I -[u] × sin ∆θ + [u] 2 × (1 -cos ∆θ) , (347)
which is a closed-form solution.

This solution corresponds to a rotation matrix, Φ = R{-u∆θ} = R{ω∆t} , according to the Rodrigues rotation formula (77), a result that could be obtained by direct inspection of (345) and recalling (69). Let us therefore write this as the final closed-form result, Φ = R{ω∆t} .

(348)

B.2 Simplified IMU example

Consider the simplified, IMU driven system with error-state dynamics governed by,

δ p = δv (349a) δ v = -R [a] × δθ (349b) δ θ = -[ω] × δθ , (349c)
where (a, ω) are the IMU readings, and we have obviated gravity and sensor biases. This system is defined by the state vector and the dynamic matrix,

x =   δp δv δθ   A =   0 P v 0 0 0 V θ 0 0 Θ θ   . ( 350
) 0 0 Θ k θ   (355)
We can observe that the terms in the increasing powers of A have a fixed part and an increasing power of Θ θ . These powers can lead to closed form solutions, as in the previous section.

Let us partition the matrix Φ as follows,

Φ =   I Φ pv Φ pθ 0 I Φ vθ 0 0 Φ θθ   , (356)
and let us advance step by step, exploring all the non-zero blocks of Φ one by one.

Trivial diagonal terms Starting by the two upper terms in the diagonal, they are the identity as shown.

Rotational diagonal term Next is the rotational diagonal term Φ θθ , relating the new angular error to the old angular error. Writing the full Taylor series for this term leads to

Φ θθ = ∞ k=0 1 k! Θ k θ ∆t k = ∞ k=0 1 k! [-ω] k × ∆t k , (357)
which corresponds, as we have seen in the previous section, to our well-known rotation matrix,

Φ θθ = R{ω∆t} . (358)

B.3 Full IMU example

In order to give means to generalize the methods exposed in the simplified IMU example, we need to examine the full IMU case from a little closer. Consider the full IMU system (237), which can be posed as

δ x = Aδx + Bw , (375)
whose discrete-time integration requires the transition matrix

Φ = ∞ k=0 1 k! A k ∆t k = I + A∆t + 1 2 A 2 ∆t 2 + . . . , (376)
which we wish to compute. The dynamic matrix A is block-sparse, and its blocks can be easily determined by examining the original equations ( 237),

A =         0 P v 0 0 0 0 0 0 V θ V a 0 V g 0 0 Θ θ 0 Θ ω 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         . (377)
As we did before, let us write a few powers of A,

A 2 =         0 0 P v V θ P v V a 0 P v V g 0 0 V θ Θ θ 0 V θ Θ ω 0 0 0 Θ 2 θ 0 Θ θ Θ ω 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         A 3 =         0 0 P v V θ Θ θ 0 P v V θ Θ ω 0 0 0 V θ Θ 2 θ 0 V θ Θ θ Θ ω 0 0 0 Θ 3 θ 0 Θ 2 θ Θ ω 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         A 4 =         0 0 P v V θ Θ 2 θ 0 P v V θ Θ θ Θ ω 0 0 0 V θ Θ 3 θ 0 V θ Θ 2 θ Θ ω 0 0 0 Θ 4 θ 0 Θ 3 θ Θ ω 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0        
.

Basically, we observe the following,

• The only term in the diagonal of A, the rotational term Θ θ , propagates right and up in the sequence of powers A k . All terms not affected by this propagation vanish. This propagation afects the structure of the sequence {A k } in the three following aspects:

• The sparsity of the powers of A is stabilized after the 3rd power. That is to say, there are no more non-zero blocks appearing or vanishing for powers of A higher than 3.

• The upper-left 3×3 block, corresponding to the simplified IMU model in the previous example, has not changed with respect to that example. Therefore, its closed-form solution developed before holds.

• The terms related to the gyrometer bias error (those of the fifth column) introduce a similar series of powers of Θ θ , which can be solved with the same techniques we used in the simplified example.

We are interested at this point in finding a generalized method to board the construction of the closed-form elements of the transition matrix Φ. Let us recall the remarks we made about the series Σ 1 and Σ 2 ,

• The subindex coincides with the lacking powers of Θ θ in each of the members of the series.

• The subindex coincides with the number of terms missing at the beginning of the series.

Taking care of these properties, let us introduce the series Σ n (X, y), defined by27

Σ n (X, y) ∞ k=n 1 k! X k-n y k = ∞ k=0 1 (k + n)! X k y k+n = y n ∞ k=0 1 (k + n)! X k y k (378)
in which the sum starts at term n and the terms lack n powers of the matrix X. It follows immediately that Σ 1 and Σ 2 respond to

Σ n = Σ n (Θ θ , ∆t) , (379)
and that Σ 0 = R{ω∆t} . We can now write the transition matrix (376) as a function of these series,

Φ =         I P v ∆t P v V θ Σ 2 1 2 P v V a ∆t 2 P v V θ Σ 3 θ ω 1 2 P v V g ∆t 2 0 I V θ Σ 1 V a ∆t V θ Σ 2 θ ω V g ∆t 0 0 Σ 0 0 Σ 1 θ ω 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I         . ( 380
)
Aδx. Closed form expressions may always be of interest, but it is unclear up to which point we should be worried about high order errors and their impact on the performance of real algorithms. This remark is particularly relevant in systems where IMU integration errors are observed (and thus compensated for) at relatively high rates, such as visual-inertial or GPS-inertial fusion schemes.

In this section we devise methods for approximating the transition matrix. They start from the same assumption that the transition matrix can be expressed as a Taylor series, and then truncate these series at the most significant terms. This truncation can be done system-wise, or block-wise.

C.1 System-wise truncation C.1.1 First order truncation: the finite differences method A typical, widely used integration method for systems of the type ẋ = f (t, x) is based on the finite-differences method for the computation of the derivative, i.e.,

ẋ lim δt→0 x(t + δt) -x(t) δt ≈ x n+1 -x n ∆t . (387)
This leads immediately to

x n+1 ≈ x n + ∆t f (t n , x n ) , (388)
which is precisely the Euler method. Linearization of the function f () at the beginning of the integration interval leads to

x n+1 ≈ x n + ∆t A x n , (389a)
where A ∂f ∂x (t n , x n ) is a Jacobian matrix. This is strictly equivalent to writing the exponential solution to the linearized differential equation and truncating the series at the linear term (i.e., the following relation is identical to the previous one),

x n+1 = e A∆t x n ≈ (I + ∆t A) x n . ( 389b
)
This means that the Euler method (App. A.1), the finite-differences method, and the firstorder system-wise Taylor truncation method, are all the same. We get the approximate transition matrix, Φ ≈ I + ∆tA . (392)

C.1.2 N-th order truncation

Truncating at higher orders will increase the precision of the approximated transition matrix. A particularly interesting order of truncation is that which exploits the sparsity of the result to its maximum. In other words, the order after which no new non-zero terms appear. For the simplified IMU example of Section B.2, this order is 2, resulting in

Φ ≈ I + A∆t + 1 2 A 2 ∆t 2 =   I I∆t -1 2 R [a] × ∆t 2 0 I -R [a] × (I -1 2 [ω] × ∆t)∆t 0 0 R{ω∆t}   . (393)
In the full IMU example of Section B.3, the is order 3, resulting in

Φ ≈ I + A∆t + 1 2 A 2 ∆t 2 + 1 6 A 3 ∆t 3 , (394)
whose full form is not given here for space reasons. The reader may consult the expressions of A, A 2 and A 3 in Section B.3.

C.2 Block-wise truncation

A fairly good approximation to the closed forms previously explained results from truncating the Taylor series of each block of the transition matrix at the first significant term. That is, instead of truncating the series in full powers of A, as we have just made above, we regard each block individually. Therefore, truncation needs to be analyzed in a per-block basis. We explore it with two examples. For the simplified IMU example of Section B.2, we had series Σ 1 and Σ 2 , which we can truncate as follows

Σ 1 = I∆t+ 1 2 Θ θ ∆t 2 + • • • ≈ I∆t Σ 2 = 1 2 I∆t 2 + 1 3! Θ θ ∆t 3 + • • • ≈ 1 2 I∆t 2 . ( 395
)
(396)

This leads to the approximate transition matrix

Φ ≈   I I∆t -1 2 R [a] × ∆t 2 0 I -R [a] × ∆t 0 0 R(ω∆t)   , (397)
which is more accurate than the one in the system-wide first-order truncation above (because of the upper-right term which has now appeared), yet it remains easy to obtain and compute, especially when compared to the closed forms developed in Section B. Again, observe that we have taken the closed-form for the lowest term, i.e., Φ θθ = R(ω∆t) .

In the general case, it suffices to approximate each Σ n except Σ 0 by the first term of its series, i.e.,

Σ 0 = R{ω∆t} , Σ n>0 ≈ 1 n! I∆t n . ( 398
)
For the full IMU example, feeding the previous Σ n into (380) yields the approximated transition matrix,

Φ ≈         I I∆t -1 2 R [a] × ∆t 2 -1 2 R∆t 2 1 3! R [a] × ∆t 3 1 2 I∆t 2 0 I -R [a] × ∆t -R∆t 1 2 R [a] × ∆t 2 I∆t 0 0 R{ω∆t} 0 -I∆t 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I         (399)
with (see ( 237)) a = a ma b , ω = ω mω b , R = R{q} , and where we have substituted the matrix blocks by their appropriate values (see also (237)),

P v = I , V θ = -R [a] × , V a = -R , V g = I , Θ θ = -[ω] × , Θ ω = -I
A slight simplification of this method is to limit each block in the matrix to a certain maximum order n. For n = 1 we have,

Φ ≈         I I∆t 0 0 0 0 0 I -R [a] × ∆t -R∆t 0 I∆t 0 0 R{ω∆t} 0 -I∆t 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I         , ( 400
)
which is the Euler method, whereas for n = 2,

Φ ≈         I I∆t -1 2 R [a] × ∆t 2 -1 2 R∆t 2 0 1 2 I∆t 2 0 I -R [a] × ∆t -R∆t 1 2 R [a] × ∆t 2 I∆t 0 0 R{ω∆t} 0 -I∆t 0 0 0 0 I 0 0 0 0 0 0 I 0 0 0 0 0 0 I         . ( 401
)
For n ≥ 3 we have the full form (399).

E.1 Noise and perturbation impulses

One is oftentimes confronted (for example when reusing existing code or when interpreting other authors' documents) with EKF prediction equations of a simpler form than those that we used here, namely,

P n+1 = F x P n F x + Q . ( 440
)
This corresponds to the general discrete-time dynamic system,

δx n+1 = F x δx n + i (441)
where i ∼ N {0, Q}

is a vector of random (white, Gaussian) impulses that are directly added to the state vector at time t n+1 . The matrix Q is simply considered the impulses covariances matrix. From what we have seen, we should compute this covariances matrix as follows,

Q = ∆t 2 B U c B + ∆t C W c C . (443)
In the case where the impulses do not affect the full state, as it is often the case, the matrix Q is not full-diagonal and may contain a significant amount of zeros. One can then write the equivalent form δx n+1 = F x δx n + F i i (444)

with i ∼ N {0, Q i } , (445)
where the matrix F i simply maps each individual impulse to the part of the state vector it affects to. The associated covariance Q i is then smaller and full-diagonal. Please refer to the next section for an example. In such case the ESKF time-update becomes δx n+1 = F x δx n (446)

P n+1 = F x P n F x + F i Q i F i . (447)
Obviously, all these forms are equivalent, as it can be seen in the following double identity for the general perturbation Q,

F i Q i F i = Q = ∆t 2 B U c B + ∆t C W c C . (448)

E.2 Full IMU example

We study the construction of an error-state Kalman filter for an IMU. The error-state system is defined in (237) and involves a nominal state x, an error-state δx, a noisy control signal u m = u + ũ and a perturbation w, specified by,

x =         p v q a b ω b g         , δx =         δp δv δθ δa b δω b δg         , u m = a m ω m , ũ = ã ω , w = a w ω w (449)
In a model of an IMU like the one we are considering throughout this document, the control noise corresponds to the additive noise in the IMU measurements. The perturbations affect the biases, thus producing their random-walk behavior. The dynamic, control and perturbation matrices are (see ( 427), ( 377) and ( 237)),

A =        
0 P v 0 0 0 0 0 0 V θ V a 0 V g 0 0 Θ θ 0 Θ ω 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

        , B =         0 0 -R 0 0 -I 0 0 0 0 0 0         , C =         0 0 0 0 0 0 I 0 0 I 0 0         (450)
In the regular case of IMUs with accelerometer and gyrometer triplets of the same kind on the three axes, noise and perturbations are isotropic. Their standard deviations are specified as scalars as follows

The system evolves with sampled measures at intervals ∆t, following (434-437), where the transition matrix F x = Φ can be computed in a number of ways -see previous appendices.

E.2.1 Noise and perturbation impulses

In the case of a perturbation specification in the form of impulses i, we can re-define our system as follows,

δx n+1 = F x (x n , u m )•δx n + F i •i (453)
with the nominal-state, error-state, control, and impulses vectors defined by,

x =         p v q a b ω b g         , δx =         δp δv δθ δa b δω b δg         , u m = a m ω m , i =     v i θ i a i ω i     , (454)
90 the transition and perturbations matrices defined by,

F x = Φ = e A∆t , F i =        
0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I 0 0 0 0

        , (455)
and the impulses variances specified by

i ∼ N {0, Q i } , Q i =     σ 2 ã∆t 2 I 0 σ 2 ω∆t 2 I σ 2 aw ∆tI 0 σ 2 ωw ∆tI     . (456)
The trivial specification of F i may appear surprising given especially that of B in (450). What happens is that the errors are defined isotropic in Q i , and therefore -Rσ 2 I(-R) = σ 2 I and -Iσ 2 I(-I) = σ 2 I, leading to the expression given for F i . This is not possible when considering non-isotropic IMUs, where a proper Jacobian F i = B C should be used together with a proper specification of Q i .

We can of course use full-state perturbation impulses,

δx n+1 = F x (x n , u m )•δx n + i (457) with i =         0 v i θ i a i ω i 0         , i ∼ N {0, Q} , Q =         0 σ 2 ã∆t 2 I 0 σ 2 ω∆t 2 I σ 2 aw ∆tI 0 σ 2 ωw ∆tI 0         . (458)
Bye bye.

  Figure 1: Rotation of a vector x, by an angle φ, around the axis u. See text for details.

Figure 2 :

 2 Figure 2: Exponential maps of the rotation matrix.

Figure 3 :

 3 Figure 3: Exponential maps of the quaternion.

Figure 6 :

 6 Figure 6: Quaternion interpolation in the unit sphere of R 4 , and a frontal view of the situation on the rotation plane π of R 4 .

Figure 8 :

 8 Figure8: Rotation in R 3 . A rotation of a vector x around an axis u describes a circumference in a plane orthogonal to the axis. The component of x parallel to the axis does not move, and is represented by the small red dot on the axis. The sketch on the right illustrates the radically different behaviors of the rotating point on the plane and axis subspaces.

Figure 10 :

 10 Figure 10: Quaternion rotation in R 4 . Two chained isoclinic rotations, one left (with equal half-angles), and one right (with opposing half-angles), produce a pure rotation by the full angle in only one of the invariant planes.

Figure 11 :

 11 Figure11: The right Jacobian J r = ∂δφ/∂δθ maps variations δθ around the parameter θ into variations δφ on the vector space tangent to the manifold at the point Exp θ.

Figure 13 :

 13 Figure 13: Integration schemes for two consecutive time steps (gray and black arrow sets), where variables sharing the same time stamp have been organized in columns. Left: forward integration. Center: midward and first-order integrations. Right: backward integration.

  -R [a] × ∆t 0 0 I -[ω∆t] ×   .(391)However, we already know from Section B.1 that the rotational term has a compact, closedform solution, Φ θθ = R(ω∆t) . It is convenient to re-write the transition matrix according to

  σ ã [m/s 2 ] , σ ω [rad/s] , σ aw [m/s 2 √ s] , σ ωw [rad/s

Table 1 :

 1 The rotation matrix and the quaternion for representing SO(3).

Table 2 :

 2 Hamilton vs. JPL quaternion conventions with respect to the 4 binary choices

Table 3 :

 3 All variables in the error-state Kalman filter.

	Magnitude	True Nominal Error	Composition	Measured Noise
	Full state ( 1 )	x t	x	δx	x t = x ⊕ δx
	Position	p t	p	δp	p t = p + δp

Table 4 :

 4 Algorithm modifications related to the definition of the orientation errors.

	Context	Item	local angular error global angular error
	Error composition	q

Our choice for the (w, x, y, z) subscripts notation comes from the fact that we are interested in the geometric properties of the quaternion in the 3D Cartesian space. Other texts often use alternative subscripts such as (0, 1,

2, 3) or (1, i, j, k), perhaps better suited for mathematical interpretations.

The skew-operator can be found in the literature in a number of different names and notations, either related to the cross operator ×, or to the 'hat' operator ∧ , so that all the forms below are equivalent,[a] × ≡ [a × ] ≡ a× ≡ a × ≡ [a] ≡ a ≡ a ∧ .

We remind that cos θ = 1θ 2 /2! + θ

/4! -• • • , and sin θ = θθ 3 /3! + θ

/5! -• • • .

We remind that arctan x = xx 3 /3 + x

/5 -• • • , and arctan(y, x) ≡ arctan(y/x).

The effect of the double cover needs to be considered when performing interpolation in the space of rotations. This is however easy, as we will see in Section 2.7.

Let Q 1 and Q 2 be orthogonal, and buildQ = Q 1 Q 2 . Then Q Q = Q 2 Q 1 Q 1 Q 2 = Q 2 IQ 2 = I.

Notice that reflections satisfy |R| = det(R) = -1, and do not form a group since |R 1 R 2 | = 1 = -1.

See footnote 6 for O(3) and add this for SO(3):let |R 1 | = |R 2 | = 1, then |R 1 R 2 | = |R 1 | |R 2 | = 1.

Put otherwise: a rotation cannot become a reflection through a continuous transformation.

The obvious notation ambiguity between the exponential maps R = Exp(v) and q = Exp(v) is easily resolved by the context: at occasions it is just the type of the returned value, R or q; other times it is the presence or absence of the quaternion product ⊗.

We can use here either the maps log() and exp(), or their capitalized forms Log() and Exp(). The involved factor 2 in the resulting angles is finally irrelevant as it cancels out in the final formula.

The angle ∆θ = arccos(q 0 q 1 ) is the angle between the two quaternion vectors in Euclidean 4-space, not the real rotated angle in 3D space, which from (125) is ∆φ = Log(q * 0 ⊗ q 1 ) . See Section 2.4.6 for further details.

Given the two invariant planes, we arbitrarily select their orientations so that we can associate positive and negative rotation angles in them.

This can be checked by extracting the eigenvalues of the isoclinic rotation matrices: they are formed by pairs of conjugate complex numbers with a phase equal to θ/2.

Let it suffice to say that q ⊗ x ⊗ q * works for rotations if q is a unit quaternion. In fact, the product q v ⊗ x ⊗ q v produces reflections (not rotations!) in R 3 if q v is a unit pure quaternion. Finally, the product q ⊗ x ⊗ q, with q a unit non-pure quaternion, exhibits no remarkable properties.

Other common denominations for the {global, local} frames are {parent, child} and {world, body}. The first one is convenient when more than two frames are involved in a system (e.g. the frames of each moving link in a humanoid robot); the second one is convenient for a solid vehicle body (e.g. a plane, a car) moving in a unique reference frame identified as the world.

Notice also from (225) that this term would always vanish if the quaternion product were commutative, which is not.

It is common practice to neglect the Earth's rotation rate ω E in the rotational kinematics described in (231), which would otherwise be ω m = ω t + R t ω E + ω bt + ω n . Considering a non-null Earth rotation rate is, in the vast majority of practical cases, unjustifiably complicated. However, we notice that when employing high-end IMU sensors with very small noises and biases, a value of ω E = 15 • /h ≈ 7.3•10 -5 rad/s might become directly measurable; in such cases, in order to keep the IMU error model valid, the rate ω E should not be neglected in the formulation.

This assumption cannot be made in cases where the three XY Z accelerometers are not identical.

x ∼ N {µ, Σ} means that x is a Gaussian random variable with mean and covariances matrix specified by µ and Σ.

We give the simplest form of the covariance update, P ← (I-KH)P. This form is known to have poor numerical stability, as its outcome is not guaranteed to be symmetric nor positive definite. The reader is free to use more stable forms such as 1) the symmetric form P ← P -K(HPH + V)K and 2) the symmetric and positive Joseph form P ← (I -KH)P(I -KH) + KVK .

Note that, being X a square matrix that is not necessarily invertible (as it is the case for X = Θ θ ), we are not allowed to rearrange the definition of Σ n with Σ n = X -n ∞ k=n 1 k! (yX) k .

The true-state kinematics

The true kinematic equations are

Here, the true acceleration a t and angular rate ω t are obtained from an IMU in the form of noisy sensor readings a m and ω m in body frame, namely 23

with R t R{q t }. With this, the true values can be isolated (this means that we have inverted the measurement equations),

Substituting above yields the kinematic system

Eliminating the second order terms and reorganizing some cross-products (with [a]

then, recalling ( 240) and ( 241),

which after proper rearranging leads to the dynamics of the linear velocity error,

To further clean up this expression, we can often times assume that the accelerometer noise is white, uncorrelated and isotropic 24 ,

that is, the covariance ellipsoid is a sphere centered at the origin, which means that its mean and covariances matrix are invariant upon rotations (Proof:

Then we can redefine the accelerometer noise vector, with absolutely no consequences, according to

Equation (237c): The orientation error. We wish to determine δ θ, the dynamics of the angular errors. We start with the following relations

which are the true-and nominal-definitions of the quaternion derivatives.

As we did with the acceleration, we group large-and small-signal terms in the angular rate for clarity,

7.1 System kinematics in continuous time 7.1.1 The true-and nominal-state kinematics

True and nominal kinematics do not involve errors and their equations are unchanged.

The error-state kinematics

We start by writing the equations of the error-state kinematics, and proceed afterwards with comments and proofs.

where, again, all equations except those of δ v and δ θ are trivial. The non-trivial expressions are developed below.

Equation (294b):

The linear velocity error. We wish to determine δ v, the dynamics of the velocity errors. We start with the following relations

where ( 295) is the small-signal approximation of R t using a globally defined error, and in (296) we introduced a B and δa B as the large-and small-signal accelerations in body frame, defined in( 240) and ( 241), as we did for the locally-defined case. We proceed by writing the expression (234b) of vt in two different forms (left and right developments), where the terms O( δθ 2 ) have been ignored,

This leads after removing Ra B + g from left and right to

Eliminating the second order terms and reorganizing some cross-products (with

where

that is, the number of iterations (the order of the method) is s, the averaging weights are defined by b i , the evaluation time instants by c i , and the slopes k i are determined using the values a ij . Depending on the structure of the terms a ij , one can have explicit or implicit RK methods.

• In explicit methods, all k i are computed sequentially, i.e., using only previously computed values. This implies that the matrix [a ij ] is lower triangular with zero diagonal entries (i.e., a ij = 0 for j ≥ i). Euler, midpoint and RK4 methods are explicit.

• Implicit methods have a full [a ij ] matrix and require the solution of a linear set of equations to determine all k i . They are therefore costlier to compute, but they are able to improve on accuracy and stability with respect to explicit methods.

Please refer to specialized documentation for more detailed information.

B Closed-form integration methods

In many cases it is possible to arrive to a closed-form expression for the integration step.

We consider now the case of a first-order linear differential equation,

that is, the relation is linear and constant over the interval. In such cases, the integration over the interval [t n , t n + ∆t] results in

where Φ is known as the transition matrix. The Taylor expansion of this transition matrix is

When writing this series for known instances of A, it is sometimes possible to identify known series in the result. This allows writing the resulting integration in closed form. A few examples follow.

with

Its integration with a step time ∆t is

The transition matrix Φ admits a Taylor development (343), in increasing powers of A∆t. We can write a few powers of A to get an illustration of their general form,

be solved by adding and substracting the missing terms, and substituting the full series by its closed form. We obtain

which is a closed-form solution valid if ω = 0. Therefore we can finally write

Position-vs-angle term Let us finally board the term Φ pθ . Its Taylor series is,

We factor out the constant terms and get,

where we note the subindex "2" in Σ 2 admits the following interpretation:

• Two powers of Θ θ are missing in each term of the series,

• The first two terms of the series are missing.

Again, we use (365) to increase the exponents of Θ θ , yielding

We substitute the incomplete series by its closed form,

which leads to the final result

Our problem has now derived to the problem of finding a general, closed-form expression for Σ n . Let us observe the closed-form results we have obtained so far,

In order to develop Σ 3 , we need to apply the identity ( 365) twice (because we lack three powers, and each application of (365) increases this number by only two), getting

which leads to

By careful inspection of the series Σ 0 . . . Σ 3 , we can now derive a general, closed-form expression for Σ n , as follows,

The final result for the transition matrix Φ follows immediately by substituting the appropriate values of Σ n , n ∈ {0, 1, 2, 3}, in the corresponding positions of (380).

It might be worth noticing that the series now appearing in these new expressions of Σ n have a finite number of terms, and thus that they can be effectively computed. That is to say, the expression of Σ n is a closed form as long as n < ∞, which is always the case. For the current example, we have n ≤ 3 as can be observed in (380).

C Approximate methods using truncated series

In the previous section, we have devised closed-form expressions for the transition matrix of complex, IMU-driven dynamic systems written in their linearized, error-state form δ x =

D The transition matrix via Runge-Kutta integration

Still another way to approximate the transition matrix is to use Runge-Kutta integration. This might be necessary in cases where the dynamic matrix A cannot be considered constant along the integration interval, i.e., ẋ(t) = A(t)x(t) .

(402)

Let us rewrite the following two relations defining the same system in continuous-and discrete-time. They involve the dynamic matrix A and the transition matrix Φ,

These equations allow us to develop ẋ(t n + τ ) in two ways as follows (left and right developments, please note the tiny dots indicating the time-derivatives),

Here, (405) comes from noticing that ẋ(t n ) = ẋn = 0, because it is a sampled value. Then,

which is the same ODE as ( 403), now applied to the transition matrix instead of the state vector. Mind that, because of the identity x(t n ) = Φ tn|tn x(t n ), the transition matrix at the beginning of the interval, t = t n , is always the identity,

with

D.1 Error-state example

Let us consider the error-state Kalman filter for the non-linear, time-varying system

where x t denotes the true state, and u is a control input. This true state is a composition, denoted by ⊕, of a nominal state x and the error state δx,

where the error-state dynamics admits a linear form which is time-varying depending on the nominal state x and the control u, i.e.,

that is, the error-state dynamic matrix in (402) has the form A(t) = A(x(t), u(t)). The dynamics of the error-state transition matrix can be written,

In order to RK-integrate this equation, we need the values of x(t) and u(t) at the RK evaluation points, which for RK4 are {t n , t n + ∆t/2, t n + ∆t}. Starting by the easy ones, the control inputs u(t) at the evaluation points can be obtained by linear interpolation of the current and last measurements,

The nominal state dynamics should be integrated using the best integration practicable. For example, using RK4 integration we have,

which gives us the estimates at the evaluation points,

We notice here that x(t n + ∆t/2) = xn+x n+1 2 , the same linear interpolation we used for the control. This should not be surprising given the linear nature of the RK update.

Whichever the way we obtained the nominal state values, we can now compute the RK4 matrices for the integration of the transition matrix,

which finally lead to,

E Integration of random noise and perturbations

We aim now at giving appropriate methods for the integration of random variables within dynamic systems. Of course, we cannot integrate unknown random values, but we can integrate their variances and covariances for the sake of uncertainty propagation. This is needed in order to establish the covariances matrices in estimators for systems that are of continuous nature (and specified in continuous time) but estimated in a discrete manner.

Consider the continuous-time dynamic system,

where x is the state vector, u is a vector of control signals containing noise ũ, so that the control measurements are u m = u + ũ, and w is a vector of random perturbations. Both noise and perturbations are assumed white Gaussian processes, specified by,

where the super-index • c indicates a continuous-time uncertainty specification, which we want to integrate.

There exists an important difference between the natures of the noise levels in the control signals, ũ, and the random perturbations, w:

• On discretization, the control signals are sampled at the time instants n∆t, having u m,n u m (n∆t) = u(n∆t) + ũ(n∆t). The measured part is obviously considered constant over the integration interval, i.e., u m (t) = u m,n , and therefore the noise level at the sampling time n∆t is also held constant, ũ(t) = ũ(n∆t) = ũn , n∆t < t < (n + 1)∆t .

(426)

• The perturbations w are never sampled.

As a consequence, the integration over ∆t of these two stochastic processes differs. Let us examine it.

The continuous-time error-state dynamics (424) can be linearized to

with

and integrated over the sampling period ∆t, giving,

which has three terms of very different nature. They can be integrated as follows:

1. From App. B we know that the dynamic part is integrated giving the transition matrix,

where Φ = e A∆t can be computed in closed-form or approximated at different levels of accuracy.

2. From (426) we have

which means that the measurement noise, once sampled, is integrated in a deterministic manner because its behavior inside the integration interval is known.

3. From Probability Theory we know that the integration of continuous white Gaussian noise over a period ∆t produces a discrete white Gaussian impulse w n described by

We obsereve that, contrary to the measurement noise just above, the perturbation does not have a deterministic behavior inside the integration interval, and hence it must be integrated stochastically. It is important and illustrative here to observe the different effects of the integration interval, ∆t, on the three terms of the covariance update (439): the dynamic error term is exponential, the measurement error term is quadratic, and the perturbation error term is linear.