e-mail comments/suggestions
Quaternions are primarily used as a representation of 3D rotations. They are convenient to use due to their compact size and their nice geometric properties. Most operations on rotations can be expressed by means of quaternion representation, which generally offers a more intuitive point of view.
This section presents a quick construction of the quaternions together with their basic properties (product, inverse, norm, etc…).
Let us start by constructing the complex numbers as real \(\small{2 \times 2}\) matrices of the form:
\[\mat{ x & -y \\ y & x \\ }\]The above set is a \(2\)-dimensional vector space, whose basis is denoted by \((1, i)\). Under this convention, one can readily check that \(i^2 = -1\) as expected.
This space is closed under matrix multiplication, and one can show that this space, together with this (commutative) multiplication, form a field, the complex numbers, denoted by \(\CC\). The complex conjugate \(\bar{x + iy} = x - i y\) corresponds to the matrix transposition:
\[\bar{x + iy} \simeq \mat{ x & -y \\ y & x \\ }^T\]Let us now repeat this construction, this time using complex matrices of the form:
\[\mat{c & d \\ -\bar{d} & \bar{c}}\](TODO: check the order) where \(c = w + i x\) and \(d = y + iz\). We end up with a 2-dimensional complex vector-space, that is also a 4-dimensional real vector space, given by real \(4 \times 4\) matrices of the form:
\[\mat{ w & -x & y & -z \\ x & w & z & y \\ -y & -z & w & x \\ z & -y & -x & w \\ }\]in which we recognize \(2\times 2\) complex blocks. This space is again closed under matrix multiplication, but this time multiplication is no longer commutative. Together with the matrix product, this space forms the quaternion group, denoted by \(\HH\).
We will generally use the real vector space structure, and identify a quaternion with its coordinates \((w, x, y, z)\) in the canonical basis.
The 4 basis vectors for the real vector space structure are commonly denoted by \((1, i, j, k)\) and satisfy the following formula:
\[i^2 = j^2 = k^2 = ijk = -1\]This formula is actually sufficient to define the quaternions, and was carved on a bridge in Dublin by Hamilton in October 1843 (cf Wikipedia). In the remaining of these notes, we will be using the following notation:
\[q = (w, x, y, z) =: w + v\]where \(v\) is an imaginary quaternion with coordinates \((0, x, y, z)\) that we will happily identify with the corresponding vector in \(\RR^3\) when needed. \(w\) is called the real part, and \(v\) the imaginary part, just like for complex numbers.
The quaternion product follows from the matrix representation product. In coordinates, it has the following expression:
\[ab = (w_a + v_a) (w_b + v_b) = w_a w_b - v_a^T v_b \ \ + \ \ w_a v_b + v_a w_b + v_a\times v_b\]In particular, for pure imaginary quaternions \(x\) and \(y\) this gives:
\[x y = x \times y - x^Ty\]This is the formula you’ll want to remember, since almost everything else follows from it.
The left multiplication by \(a\) can be put in matrix form as:
\[a b = \underbrace{\mat{w_a & -v_a^T\\ v_a & w_a I + \hat{v}_a}}_{L_a} \mat{w_b\\v_b}\]Similarly, the right multiplication by \(b\) can be put in matrix form as:
\[a b = \underbrace{\mat{w_b & -v_b^T\\ v_b & w_b I - \hat{v}_b}}_{R_b} \mat{w_a\\v_a}\]The associativity of the quaternion product implies that \(L_a\) and \(R_b\) commute for any quaternions \(a, b\). Finally, one can easily check the following identities:
\[\begin{align} R_{a + b} &= R_a + R_b \\ L_{a + b} &= L_a + L_b \\ R_{\lambda q} &= \lambda R_q\\ L_{\lambda q} &= \lambda L_q\\ \end{align}\]As for the complex numbers, the conjugate corresponds to the transpose in the matrix representation. This immediately implies that:
\[\bar{p q} = \bar{q} \ \bar{p}\]As for complex numbers, \(\bar{1}=1\), \(\bar{i} = -i\) and furthermore \(\bar{j} = -j\), \(\bar{k} = -k\) so that conjugation of an arbirary quaternion is simply given by:
\[\bar{q} = w_q - v_q\]The real and imaginary parts of a quaternion may be obtained using conjugation:
\[w_q = \frac{q + \bar{q}}{2}\] \[v_q = \frac{q - \bar{q}}{2}\]In matrix form, the conjugation operator is given by:
\[H = \matrix{1 & 0 \\ 0 & -I}\]One can easily check that the left/right translation matrices satisfy the following identities:
\[\begin{align} L_{\bar{q}} &= L_q^T\\ R_{\bar{q}} &= R_q^T\\ \end{align}\]As for the complex numbers, conjugation induces a norm over the quaternions (through the Frobenius norm on the matrix representation), which coincides with the Euclidean norm on \(\RR^4\):
\[|q|^2 = q \bar{q} = \bar{q} q = ||q||_{\RR^4}^2 \geq 0\]As for the complex numbers, the norm is multiplicative:
\[|ab| = |a|\ |b|\]From the norm definition, a simple expression can be obtained for the quaternion inverse:
\[\inv{q} = \frac{\bar{q}}{|q|^2}\]The unit sphere \(S^3\) is closed under quaternion multiplication and inverse, and it contains the identity element \(1\): this makes it a multiplicative subgroup of \(\HH\). As we shall see, this space is particularly interesting for representing rotations. The unit real 4-sphere with quaternion multiplication is isomorphic to the quaternion representation using unit complex coefficients, i.e. the group \(SU(2)\).
As the inverse image of \(0\) by the smooth function \(f(q) = \norm{q}^2 - 1\), the unit hypersphere \(S^3\) is a closed, \(3\)-dimensional smooth sub-manifold of \(\RR^4\). It is also compact and simply connected, meaning that every smooth curve from \(S^1\) to \(S^3\) can be shrunk into a point.
\(S^3\) is a multiplicative subgroup of \(\HH\) and since both multiplication and the inverse are smooth, it also forms a Lie group. The Lie group structure provides the connection with rotations in \(\RR^3\) through the adjoint representation.
The tangent space at the identity \(1\) is the space of imaginary quaternions, which we identify with \(\RR^3\). The inner automorphism \(\Psi_q\):
\[\Psi_q: S^3 \to S^3 \quad h \mapsto q h \bar{q}\]is sometimes known as the conjugation by \(q\), and its derivative at \(1\) is the adjoint of the Lie group:
\[\Ad_q = \dd \Psi_q(1): \RR^3 \to \RR^3 \quad x \mapsto qx\bar{q}\]From the multiplicative property of the norm, we see that \(\Ad_g\) is an isometry of \(\RR^3\), so it is either a rotation or a reflection (or a mix of the two). Since \(S^3\) is connected and \(1 \in S^3\) (and \(\Ad_1 = I_3\) which is a rotation) \(\Ad_g\) has to be a rotation: there is a smooth path from \(1\) to any quaternion so that \(\det\block{\Ad_g}\) has to remain \(1\).
\(\Ad\) provides a group homomorphism between \(S^3\) and \(SO(3)\) (the adjoint representation):
\[\Ad: S^3 \to SO(3)\]so that multiplying quaternions corresponds to composing rotations. However, it is not injective (the representation is not faithful): both \(q\) and \(-q\) share the same adjoint. This \(2\)-to-\(1\) relationship is known as the double covering of \(SO(3)\) by \(S^3\).
The derivative of the adjoint map at the identity provides the Lie algebra structure on \(\RR^3\), noted \(\alg{s^3}\):
\[\ad = \dd \Ad_1: \RR^3 \mapsto \alg{so(3)}\]where the Lie bracked is given by:
\[[x, y] = \ad(x)y\]The actual derivation of \(\Ad\) at the identity gives an expression as the quaternion commutator:
\[[x, y] = xy - yx\]In coordinates, the Lie bracked is twice the cross product:
\[[x, y] = 2 x \times y\]While we’re at it, \(\ad\) also provides a Lie algebra isomorphism between \(\alg{s^3}\ (\simeq \RR^3)\) and \(\alg{so(3)}\), the space of \(3 \times 3\) skew-symmetric matrices:
\[\ad(x) = 2 \hat{x}\]where \(\hat{\omega} = \mat{0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0}\) is the cross-product anti-symmetric matrix, such that \(\hat{\omega}x = \omega \times x\). The Killing form for \(\alg{s^3}\) satisfies:
\[\left\langle \hat{\omega}_1, \hat{\omega}_2 \right\rangle = -\half \tr\block{\hat{\omega}^T_1\hat{\omega}_2} = \omega_1^T \omega_2\]Finally, as a Lie algebra isomorphism, \(\ad\) preserves the Lie bracket:
\[\ad( [x, y] ) = [\ad(x), \ad(y)]\]The classical characterization of the exponential (due to Euler) conveys pretty well what the exponential does:
\[\exp(v) = \lim_{n \to \infty} \block{ 1 + \frac{v}{n} }^n\]That is: we cut our vector \(v\) in infinitesimally small pieces to obtain elements of the group infinitesimally close to \(1\) in the direction \(v\), and multiply (compose) these elements together. As in the matrix case, one can verify that the exponential is also given by the usual Taylor series, this time using quaternion product:
\[\exp(v) = \sum_{i=0}^{\infty} \frac{v^i}{i!}\]When \(v \neq 0\), there exists a unit vector \(n\in S^2\) such that \(v = \norm{v} n\). This unit vector satisfies \(n^2 = -1\), so the above sum can be rewritten as:
\[\begin{align} \exp(v) &= \sum_i \frac{v^{2i}}{2i!} + \frac{v^{2i+1}}{2i+1!} \\ &= \sum_i (-1)^i \frac{\norm{v}^{2i}}{2i!} + (-1)^{i} \frac{\norm{v}^{2i+1}}{2i+1!} n \\ &= \cos \norm{v} + \sin \norm{v} n \end{align}\]The injectivity radius is \(2\pi\), but we start getting the same rotations again after a radius of \(\pi\).
The logarithm is the inverse operation, defined inside the injectivity radius:
\[\log(q) = \theta n\]where:
\[\begin{align} \theta &= \arccos(w) \\ n &= \frac{v}{\norm{v}} \\ \end{align}\]Since both \(q\) and \(-q\) represent the same rotation, it is common to flip the quaternion to the positive real hemisphere prior to computing the logarithm:
\[q^+ = \begin{cases} q & \text{if} & w_q \geq 0\\ -q & \text{if} & w_q < 0\\ \end{cases}\]Doing so ensures the arc-cosine is well defined, and that the rotation interpolated along the logarithm by \(\exp\block{\alpha \log(q)},\,\alpha \in [0, 1]\) takes the short way.
At this point, we still do not know which rotation corresponds to the conjugation by a unit quaternion. From the Taylor series, we see that \(\exp(v)\) and \(v\) commute, so that:
\[\Ad_{\exp(v)} = \exp(v)v\exp(-v) = \exp(v)\exp(-v)v = v\]and the rotation axis of \(\Ad_{\exp(v)}\) is \(v\). The one-parameter (additive) subgroup generated by \(\exp(\alpha v)\) corresponds to a one-parameter subgroup of rotations, i.e. along the same axis. Hence the rotation angle is propotional to \(\norm{v}\). Finally, \(\exp(\pi n) = -1\) for a unit vector \(n\), which corresponds to the identity rotation (or \(2\pi\)), so the rotation angle for \(\exp(v)\) is \(2\norm{v}\). Equivalently:
\[R_{n, \theta} \simeq \cos \block{\frac{\theta}{2}} + \sin\block{\frac{\theta}{2}} n\]Now hopefully that double covering story starts making sense: a half-hemisphere of \(S^3\) covers the whole \(SO(3)\), and the whole \(S^3\) covers \(SO(3)\) twice. The exponential map for \(SO(3)\) and \(S^3\) are related by:
\[\exp_{SO(3)}\block{\hat{x}} = \Ad\block{ \exp_{S^3}\block{\frac{x}{2}}}\]The derivatives (here in body-fixed coordinates) are related by:
\[\begin{align} \db \exp_{SO(3)}\block{\hat{x}}.\dd \hat{x} &= \half \ad_{\alg{s^3}}\block{\db \exp_{S^3}\block{\frac{x}{2}}}.\dd x\\ &= \block{\db \exp_{S^3}\block{\frac{x}{2}}.\dd x}^{\hat{}} \end{align}\]The standard metric in \(\RR^4\) induces a Riemannian manifold structure on \(S^3\), which plays nicely with the Lie group structure as we shall see.
The ambient metric on \(\RR^4\) induces a metric on each fiber (tangent space) of the tangent bundle \(TS^3\). This metric is left-invariant:
\[\inner{u, v}_{\RR^4} = \inner{qu, qv}\]It is also right invariant (and thus also \(\Ad\)-invariant):
\[\inner{u, v}_{\RR^4} = \inner{uq, vq}\]To see both of these, it is useful to come back to the matrix representation of quaternions, in which the metric is:
\[\inner{u, v} = \trace{u^Tv}\]and the \(q, q^T\) cancel each other nicely. Therefore, the Lie algebra metric extended to the whole manifold by left or right translation, is a Riemannian metric (the one induced by the ambient space), called a bi-invariant metric. An interesting consequence is that the Lie group exponential can be used to compute geodesics for this metric.
More generally, compact Lie groups (such as \(S^3\)) always have a bi-invariant metric.
Geodesics are locally shortest paths between points in the manifold. The fact that the Riemannian metric is bi-invariant implies that curves of the form:
\[q \exp(t x),\quad t\in\RR\]are geodesic curves (and all geodesics are of this form). We can obtain geodesic curves joining two unit quaternions \(a, b \in S^3\) as follows:
\[f(t) = a\exp\block{t \log\block{ \inv{a} b } }\]This is generally referred to as the spherical linear interpolation, or SLERP, and is used to interpolate the corresponding rotations.
The Riemannian metric induces a distance on the manifold, obtained by measuring the length of the shortest geodesic curve between two unit quaternions. The length of the geodesic between \(a\) and \(b\) (and thus the geodesic distance) is:
\[d(a, b) = \norm{\log\block{\inv{a}b}}\]Expanding a little bit, we have:
\[\begin{align} d(a, b) &= \arccos\block{w_{\inv{a}b}} \\ &= \arccos\block{ w_a w_b + v_a^T v_b } \\ \end{align}\]which is exactly the angle between 4-vectors \(a\) and \(b\), as one would expect from vectors on a unit sphere.
After working out the complete conjugation operation1, one can obtain the following expression:
\[\begin{align} q x \bar{q} &= x + 2 w \hat{v} x + 2 \hat{v}^2 x \\ &= \block{I + 2 \cos\block{\frac{\theta}{2}} \sin\block{\frac{\theta}{2}} \hat{n} + 2 \sin^2\block{\frac{\theta}{2}} \hat{n}^2} x \\ \end{align}\]Using the identity \(1 - 2 \sin^2\block{\frac{\theta}{2}}=\cos(\theta)\) and \(2 \cos\block{\frac{\theta}2} \sin \block{\frac{\theta}2}=\sin(\theta)\), we obtain the Euler-Rodrigues formula:
\[R = \cos(\theta) I + \sin(\theta) \hat{n} + \block{1 - \cos\block{\theta}}n n^T\]Or, emphasizing the action on stable subspaces:
\[R = nn^T + \cos(\theta) (I - nn^T) + \sin(\theta) \hat{n}\]From the above we easily obtain the formula for the trace of a rotation matrix:
\[\trace{R} = 1 + 2\cos(\theta)\]since \(\trace{nn^T} = 1\) and \(\tr{\hat{n}} = 0\).
Conversely, a rotation matrix can be converted back to one of its corresponding unit quaternion. From the above formula, the skew-symmetric part of a rotation \(R\) is:
\[\begin{align} \frac{R - R^T}{2} &= 2 w \hat{v} \\ &= 2 \cos\block{\frac{\theta}2} \sin \block{\frac{\theta}2}\hat{n} \\ &= \sin(\theta) \hat{n} \\ \end{align}\]which we can use to recover \(n\). But if \(\theta\) is close to \(\pi\), we cannot recover the axis from the skew-symmetric part alone, since the formula degenerates to:
\[R = I + 2\hat{n}^2 = I + 2\block{n n^T - I} = 2 n n^T - I\]which is symmetric. By adding \(I\) to \(R\), we get \(2n n^T\) which we can use to obtain \(n\):
\[n n^T = \mat{n_x^2 & n_x n_y & n_x n_z \\ n_y n_x & n_y^2 & n_y n_z \\ n_z n_x & n_z n_y & n_z^2 }\]The diagonal provides absolute values for \(n_x, n_y, n_z\). Since we’re dealing with half-turns, \(n\) is only defined up to sign so we may freely choose the one with positive \(n_x\). In this case, we get the sign for \(n_y\) using the sign of off-diagonal entries in \(nn^T\):
\[\mathrm{sign}\block{n_y} = \mathrm{sign}\block{n_x n_y}\]and similarly for \(n_z\).
As seen before, any unit quaternion may be written as:
\[q = \cos(\theta) + \sin(\theta) n\]where \(\theta \in S^1\) and \(n \in S^2\) (which gives the Hopf fibration). Let us now consider the derivative of this formula:
\[\dd q = -\sin(\theta) \dd \theta + \cos(\theta) n \dd \theta + \sin(\theta) \dd n\]To keep things readable, we introduce \(c = \cos(\theta)\) and \(s = \sin(\theta)\), and remark that \(n^T\dd n = 0\) since \(n \in S^2\), hence the pure quaternion product \(n\dd n = n \times \dd n\), and \(n^2 = -1\). The body-fixed derivative is:
\[\begin{align} \db q &= \inv{q} \dd q \\ &= (c - s n) (-s \dd \theta + c n \dd \theta + s \dd n) \\ &= -cs \dd \theta + c^2 n \dd \theta + cs \dd n + s^2 n \dd \theta -sc n^2 \dd \theta - s^2 n \dd n \\ &= n\ \dd \theta + cs\ \dd n - s^2 n\ \dd n\\ &= n\ \dd \theta + s \block{c - s n} \dd n\\ \end{align}\]By inverting this formula in the respective stable subspaces, one can obtain the following:
\[\begin{align} \dd \theta &= n^T \db q\\ \dd n &= \frac{1}{\sin(\theta)} \block{cI + s \hat{n}} \block{I - nn^T}\db q\\ &= \frac{1}{\tan \theta}\block{I - n n^T} \db q + \hat{n} \db q \\ \end{align}\]The polar decomposition is useful when computing the logarithm derivative. In polar coordinates, the logarithm is simply:
\[\log(\theta, n) = \theta n\]Thus, when defined, the derivative of the logarithm satisfies:
\[\dd \log(\theta, n) = \dd \theta \ n + \theta \dd n\]The body-fixed derivative may be obtained by plugging the orthogonal decomposition for \(\dd \theta, \dd n\) given \(\db q\). Which gives:
\[\db \log(q).\db q = nn^T \db q + \frac{\theta}{\tan(\theta)}\block{I - nn^T}\db q + \log(q)\times \db q\]Note that this formula can be extended by continuity at \(q=1\). See Bullo et al.2 for a similar formula in the case of \(SO(3)\), replacing \(\block{\theta, \log(q)}\) by \(\block{\theta/2, \log(q)/2}\) to account for the double covering.
Similarly, given a pure imaginary quaternion \(x\in\alg{s^3}\) one can obtain the polar coordinates of \(\exp(x)\) as:
\[(\theta, n) = \block{\norm{x}, \frac{x}{\norm{x}}}\]Therefore:
\[\begin{align} \dd \theta &= \frac{x^T}{\norm{x}} \dd x \\ \dd n &= \frac{1}{\norm{x}}\block{I - \frac{x x^T}{\norm{x}^2}} \dd x \\ \end{align}\]Plugging the above in the body-fixed derivative formula for polar coordinates gives:
\[\begin{align} \db \exp(x).\dd x = \db q &= n \dd \theta + s(c - sn) \dd n \\ &= nn^T \dd x + \frac{\sin(\theta)}{\theta}\block{cI - s \hat{n}} \block{I - nn^T} \dd x \\ \end{align}\]It is worth having a closer look at this formula. For a given \(\dd x\):
As seen above, the SLERP between unit quaternions \(a\) and \(b\) is expressed as:
\[f(t) = a\exp\block{t \log\block{ \inv{a} b } }\]Kim et al3 propose spline-like interpolation curves for unit quaternions \(\block{q}_{i \geq 0}\) as follows: the sline basis functions \(B_i(t)\) are first transformed to cumulative basis functions \(\tilde{B}_i(t) = B_i(t) - B_{i-1}(t)\), then the curve is obtained as:
\[f(t) = q_0 \prod_{i=1}^n \exp\block{ \tilde{B}_i(t) \log\block{q_{i-1}^{-1}q_i}}\]The authors report low acceleration/torque and easy to compute derivatives. This scheme is easily generalizable to other Lie groups.
The Euclidean mean minimizes the sum of distances squared to all the points considered. This property is easily generalized to Riemannian manifolds, provided the minimizer exists and is unique. Since \(S^3\) is compact, the existence is granted, but the unicity is not (think of averaging antipodal points).
Arsigny et al4, as others before5, describe a fast, gradient-descent-like iterative scheme to compute the geometric mean of \(n\) unit quaternions \((q_i)_{i \geq 1}\). The algorithm produces successive estimates of the geometric mean \(\mu_k\) as follows:
\[\begin{align} \mu_0 &= 1 \\ \mu_{k+1} &= \mu_k \exp\block{ \frac{1}{n} \sum_{i=0}^n \log\block{ \mu_k^{-1} q_i } } \\ \end{align}\]That is: the data is linearized in the tangent space at \(\mu_k\), then the Euclidean mean is computed, and the exponential maps the result back to the group. The algorithm stops when the tangent Euclidean mean is sufficiently close to zero. This algorithm generalizes easily to other Lie groups, but the existence might be harder to establish.
Right translations can of course be used instead, depending on the problem.
The Cayley map for imaginary quaternions has an interesting geometric interpretation:
\[\begin{align} Cay(x) &= \frac{1 - x}{1 + x} \\ &= \frac{ \block{1 - x}^2}{ \block{1 - x}\block{1 + x} } \\ &= \block{ \frac{1 - x}{\norm{1 - x}} }^2 \\ \end{align}\]That is:
Now, since \(Cay\) is its own inverse, one can interpret its computation for a unit quaternion as the following (reversed) steps:
So the square root can be expressed in terms of \(Cay\) and other simple operations. In particular, we have:
\[q^{\half} = \pi_{S^3} \block{ 1 - \frac{1 - q}{1 + q} }\]A much simpler formula can also be obtained:
\[1 - \frac{1 - q}{1 + q} = \frac{2q}{1+q} = 2 \frac{q + \norm{q}^2}{ \norm{1 + q}^2 }\]so that:
\[q^\half = \pi_{S^3} \block{q + 1}\]That is: add \(1\) to the real part, then normalize. This is a direct extension of a similar formula for computing the square root of unit complex numbers. For non-unit quaternions, one can always take the square-root of the norm separately and reduce to the case above:
\[q^\half = \sqrt{\norm{q}} \pi_{S^3} \block{q + \norm{q}}\]The quaternion product for imaginary quaternions is:
\[x y = x \times y - x^T y = \norm{x} \norm{y} \block{-\cos \theta + n \sin \theta}\]where \(\theta\) is the angle between \(x\) and \(y\). This expression, once normalized, gives us twice the rotation from \(y\) to \(-x\), so we want to consider the product \(-yx\) instead. We have just seen how to obtain square roots efficiently, so the following does what we need:
\[q_{x \mapsto y} = \pi_{S^3} \block{1 + \pi_{S^3} \block{-y x}}\]Or even simpler:
\[q_{x \mapsto y} = \pi_{S^3} \block{ \norm{y x} - y x }\]Alternatively, one can simply normalize \(-(x+y)x\) (resp. \(-y(x + y)\)), giving twice the rotation from \(x\) to \(x+y\) (resp. from \(x + y\) to \(y\)).
There is a closed-form solution to the problem of the geodesic projection of a unit quaternion to a geodesic. This problem arises in non-Euclidean generalizations6 of the Principal Component Analysis. Assuming we are given a geodesic curve starting at \(1\), with (normalized) direction \(n \in S^2\), and a unit quaternion \(q\), we consider the following optimization problem:
\[\argmin{\alpha \in \RR}\quad f(\alpha) = \half d\block{\exp\block{\alpha n}^{-1}, q}^2\]We may rewrite the objective function a little bit:
\[\begin{align} f(\alpha) &= \half \norm{ \log \block{ \exp\block{-\alpha n} q } }^2 \\ \end{align}\]We introduce \(e = \exp\block{-\alpha n}\) and \(u = e q\). The derivative, using spatial coordinates where needed, is:
\[\begin{align} \dd f(\alpha) &= -\log(u)^T \ds \log(u) \ds \exp(-\alpha n) n \\ &= -\log(u)^T n \\ \end{align}\]The first-order optimality conditions for our problem imply that:
\[\log(u)^T n = 0\]Since the logarithm will be along the imaginary part, it is equivalent to ask for the imaginary part of \(u\) to be orthogonal to \(n\):
\[v_u = w_e v_q + w_q v_e + v_e \times v_q\]\(v_e = -\sin (\alpha) n\) so we may drop the cross product which will be orthogonal to \(n\). We are left with:
\[w_e v_q + w_q v_e \quad \bot \quad n\]Or:
\[\cos(\alpha) v_q^T n - w_q \sin(\alpha) = 0\]if \(w_q \neq 0\):
\[\begin{align} \frac{v_q^T n}{w_q} &= \tan(\theta) \\ \alpha &= \theta \mod \pi \end{align}\]With \(x = \frac{v_q^T n}{w_q}\), the projected quaternion is:
\[\exp\block{\theta n} = \frac{1 + x n}{\sqrt{1 + x^2}} = \frac{1 + nn^T\frac{v_q}{w_q}}{\sqrt{1 + x^2}}\]Geometrically, this corresponds to:
This is exactly the gnomonic projection for \(S^3\), as one could expect.
if \(v_q^T n \neq 0\):
\[\begin{align} \frac{w_q}{v_q^T n} &= \tan(\theta) \\ \alpha &= \pi / 2 - \theta \mod \pi \end{align}\]TODO figure out whether Said et al7 got the same formula.
The gnomonic projection is useful for geodesic projection and is given by:
\[\pi(q) = \frac{q}{w} = \frac{q}{\cos \theta} = 1 + \tan(\theta) n\]The derivative is thus:
\[\begin{align} \dd \pi(q).\dd q &= \frac{\dd q}{w} + \frac{\sin(\theta)}{w^2}\dd \theta \\ &= n\block{1 + \tan^2(\theta)}.\dd \theta + \tan(\theta).\dd n \\ \end{align}\]where the two equivalent expressions for \(\pi\) were used on each line. Using body-fixed coordinates for \(\dd q\) and the derivative formula for polar coordinates yields in both cases:
\[\begin{align} \dd \pi(q).\db q &= \pi(q) \block{\tan(\theta) n^T \db q + \db q} \\ &= \db q + \pi(q) \times \db q + \pi(q) \pi(q)^T \db q \\ &= \block{I + \tan(\theta) \hat{n} + \tan^2(\theta) n n^T} \db q \\ \end{align}\]Here is a geometric approach to finding the Euler angles corresponding to a unit quaternion. We look for the following decomposition:
\[q \in S^3 = \exp(\theta_x e_x) \exp(\theta_y e_y) \exp(\theta_z e_z)\]for some Euler angles \(\theta_x, \theta_y, \theta_z\). The gnomonic projection \(\pi: S^3 \to \alg{s^3}\) is useful in this context since it maps geodesics to straight lines: the geodesic along \(e_z\) has body-fixed velocity \(e_z\) at \(q\), which in gnomonic projection corresponds to a line passing through \(\pi(q)\). The line direction is given by the gnomonic projection derivative formula derived above:
\[u = e_z + \pi(q)\times e_z + \pi(q) \pi(q)^T e_z\]By the same token, the geodesic along \(e_x\) corresponds to a line starting at \(0\) with direction \(e_x\), and finally the geodesic along \(e_y\) corresponds to a line starting at some point along the first geodesic projection \(p_x = \lambda_x e_x\), with a direction corresponding to a body-fixed velocity \(e_y\) at \(p_x\). This direction \(v\) is again given by the gnomonic projection derivative formula above:
\[\begin{align} v &= e_y + \underbrace{p_x \times e_y}_{\lambda_x e_x \times e_y} + p_x \underbrace{p_x^T e_y}_0 \\ &= e_y + \lambda_x e_z \\ \end{align}\]for some real \(\lambda_x\). We want the geodesic along \(e_y\) and the one along \(e_z\) to meet, so we look for some point satisfying:
\[\pi(q) + \lambda_u u = \lambda_x e_x + \lambda_v v\]Expanding the right-hand side, we obtain the following system of equations to solve for \(\lambda_x, \lambda_u, \lambda_v\):
\[\pi(q) + \lambda_u u = \lambda_x e_x + \lambda_v\block{e_y + \lambda_x e_z} = \mat{\lambda_x \\ \lambda_v \\ \lambda_x \lambda v}\]Geometrically, we are looking for the intersection of an affine line \(\pi(q) + \lambda_u u\) with the hyperbolic paraboloid \(z = x y\). Substituting \(x, y, z\) with their values from \(\pi(q) + \lambda_u u\), one finally obtains a quadratic equation in \(\lambda_u\), which should always have real solution(s). Solving for \(\lambda_u\) gives the quaternion:
\[\pi^{-1}\block{\pi(q) + \lambda_u u} = \exp(\theta_x e_x) \exp(\theta_y e_y)\]Likewise \(\lambda_x\) gives the quaternion:
\[\pi^{-1}\block{\lambda_x e_x} = \exp(\theta_x e_x)\]It is then easy to obtain the Euler angles \(\theta_x, \theta_y, \theta_z\). It is remarkable that the above method obtains the 2 intermediate quaternions using purely geometric operations. In particular, no trigonometric functions were used. A similar procedure can also be used to find Euler angles along arbitrary (orthogonal) rotation axis.
For an imaginary quaternion \(x\), one can easily check that the following formula holds:
\[\dd \block{x^n}. \dd x = n\ x^{n-1} \dd x_\parallel + \sum_{i=n-1}^{i = 0} x^{n-1} \dd x_\perp\]where \(\dd x_\parallel\) is the projection of \(\dd x\) onto \(x\) and \(\dd x_\perp = x - \dd x_\parallel\). The above is equivalent to:
\[\begin{align} \dd \block{x^{2n}}.\dd x &= (2n) x^{2n-1} \dd x_\parallel \\ \dd \block{x^{2n+1}}.\dd x &= (2n+1) x^{2n} \dd x_\parallel + x^{2n} \dd x_\perp \\ \end{align}\]Note: the quaternionic product is used. This formula can be used to compute the derivative of the exponential map:
\[\begin{align} \dd \exp(x).\dd x &= \sum_i \frac{\dd \block{x^i}}{i!} \dd x \\ &= \sum_i \frac{x^i}{i!} \dd x_\parallel + \sum_i \frac{x^{2i}}{ (2i+1)!} \dd x_\perp \\ &= \exp(x)\ \dd x_\parallel + \frac{\sin \theta}{\theta} \dd x_\perp \\ \end{align}\]Again, note that the quaternion product is used. The body-fixed derivative corresponds to the formula obtained previously.
In practice, one is generally concerned with obtaining the second derivative at the origin. The power series give:
\[\exp(x) = \sum_i \frac{x^i}{i!}\] \[\dd \exp(x).\dd x = \dd x + \frac{\dd x.x + x.\dd x}{2} + \frac{\dd x.x^2 + x.\dd x.x + x^2.\dd x}{6} + \ldots\] \[\begin{align} \dd^2 \exp(x).\dd x_2.\dd x_1 &= \frac{\dd x_1.\dd x_2 + \dd x_2.\dd x_1}{2} + O(x) \\ &= -\dd x_1^T \dd x_2 + O(x) \end{align}\]which gives at \(x = 0\):
\[\dd^2 \exp(0).\dd x_1.\dd x_2 = -\dd x_1^T \dd x_2\]As expected, the expression is symmetric and orthogonal to the Lie algebra.
To simulate rigid bodies in a stable manner, one needs to compute the second derivative of the following mapping, for a given \(y \in \RR^3\):
\[f_y: x \mapsto \Ad_{\exp(x)} y\]where \(f_y\) maps the local coordinates of a point to the absolute frame through a rotation in exponential coordinates. As before, we get:
\[\dd^2 \exp\block{x}.\hat{\dd x_2}.\hat{\dd x_1} = \frac{\hat{\dd x_1}.\hat{\dd x_2} + \hat{\dd x_2}.\hat{\dd x_1}}{2} + O\block{\hat{x}}\]From the double cross-product formula \(\hat{u} \hat{v} = vu^T - u^TvI\), we get:
\[\frac{\hat{\dd x_1}\hat{\dd x_2} + \hat{\dd x_2}\hat{\dd x_1}}{2} = \frac{\dd x_1 \dd x_2^T + \dd x_2 \dd x_1^T}{2} - \dd x_1^T\dd x_2 I\]Again, the expression is symmetric and orthogonal to the Lie algebra. Given a force \(\lambda \in \RR^3\) applied on \(f_y(0)\), the associated geometric stiffness is:
\[\begin{align} \lambda^T \dd^2 f_y (0).\dd x_1.\dd x_2 &= \lambda^T\frac{\hat{\dd x_1}.\hat{\dd x_2} + \hat{\dd x_2}.\hat{\dd x_1}}{2} y\\ &= \dd x_1^T \frac{ \hat{y} \hat{\lambda} + \hat{\lambda} \hat{y}}{2} \dd x_2 \end{align}\]Consider a reflection \(M_n\) by a plane whose normal is unit vector \(n\):
\[\begin{align} M_n(x) &= x^\parallel - x^\bot \\ &= \block{I - n n^T} x - n n^T x \\ &= \block{I - 2 n n^T} x \end{align}\]This is almost the half-turn around \(n\) as given by the Rodrigues formula:
\[\Ad_n(x) = x + 2 \hat{n}^2 x = \block{2 n n^T - I}x\]We immediately obtain:
\[\begin{align} M_n(x) &= -\Ad_n(x) \\ &= n^2 n^{-1} x n \\ &= n x n \end{align}\]Similarly, we can derive expressions for orthogonal projections:
\[x^\parallel = \frac{x - n x n}2\] \[x^\bot = \frac{x + n x n}2\]Let us consider two sets of \(n\) 3-dimensional points arranged as flat matrices \(x, y \in \mathcal{M}_{3, n}\), and look for the rotation matrix that best aligns \(x\) over \(y\):
\[\argmin{R \in SO(3)}\quad \half \norm{Rx - y}^2\]which we can rewrite as the following problem:
\[\argmin{R \in SO(3)}\quad -\tr\block{y^T R x}\]or, equivalently, as a maximization over unit quaternions:
\[\argmax{q \in S^3}\quad \sum_i \inner{y_i, q x_i \bar{q}}_{\RR^4}\]Now, we can exploit the bi-invariance of the Euclidean metric to rewrite each term of the sum as:
\[\begin{align} \inner{y_i, q x_i \bar{q}} &= \inner{y_i q, q x_i} \\ &= \inner{L_{y_i} q, R_{x_i} q} \\ &= q^T \underbrace{L_{y_i}^T R_{x_i}}_{M_i} q \\ \end{align}\]so that our initial problem is reduced to that of maximizing:
\[\argmax{q \in S^3}\quad q^T \underbrace{\block{\sum_i M_i}}_M q\]As usual, the stationary condition give:
\[\exists \lambda \in \RR, \quad Mq = \lambda q\]Since we are maximizing \(q^TMq = \lambda\), the solution quaternion is given by the largest eigenvector of matrix \(M\).
Please refer to the dedicated page.
The derivation uses the double cross-product identity: \(v \times (v \times x) = v.v^Tx - x.v^T v\) ↩
F. Bullo , R. M. Murray, Proportional Derivative (PD) Control On The Euclidean Group, European Control Conference, 1995. (Lemma 3) ↩
Kim, Myoung-Jun, Myung-Soo Kim, and Sung Yong Shin. “A general construction scheme for unit quaternion curves with simple high order derivatives.” Proceedings of the 22nd annual conference on Computer graphics and interactive techniques. ACM, 1995. ↩
Arsigny, Vincent, Xavier Pennec, and Nicholas Ayache. “Bi-invariant means in Lie groups. Application to left-invariant polyaffine transformations.” (2006). ↩
Moakher, Maher. “Means and averaging in the group of rotations.” SIAM journal on matrix analysis and applications ↩
Fletcher, P. Thomas, et al. “Principal geodesic analysis for the study of nonlinear statistics of shape.” Medical Imaging, IEEE Transactions on 23.8 (2004): 995-1005. ↩
Said, Salem, et al. “Exact principal geodesic analysis for data on so (3).” 15th European Signal Processing Conference (EUSIPCO-2007). EURASIP, 2007. ↩
e-mail comments/suggestions