Optimizing Motion Sequences with Projective Dual Quaternions

1. Introduction

Attitude control of robots, aircraft and spacecraft involves rotations about fixed gimbals, for which Euler-type decompositions are applied in order to bring the system to its desired state [1,2,3,4]. However, due to the famous gimbal lock problem [5], arising from the topological singularity in the parametrization map between the three-dimensional torus and the rotation group ${\mathbb{T}}^3\leftarrow {\mathbb{RP}}^3\cong \mathrm{SO}\left(3\right)$, as well as for sheer convenience, other representations are often preferred, such as Hamilton’s unit quaternions and Rodrigues’ vectors, which are in essence projective quaternions [6,7,8]. Similarly, for more general rigid motions, including rotational component and displacement, it is convenient to use the dual analogues of these entities and extend the usual treatment of the rotation group to the affine case by Kotelnikov’s principle of transference known from the theory of screws [9,10,11,12,13].

The problem of optimization is quite reasonable in both cases, and while for pure rotations it has been studied for a long time using both analytical and geometric tools [14,15,16,17,18,19], general rigid motion seems a bit more challenging and has been addressed mainly with machine learning algorithms [20,21,22] and enjoys less attention overall. Here we consider both in a unified geometric approach, deriving closed form solutions based on generalized Euler decompositions previously introduced in [4]. The idea is quite similar to the one introduced in [5], namely that the gimbal lock setting allows for optimizing the sum of the angles (length of the spherical path) in a rotation sequence using the redundant free parameter. This paper expands on that concept by introducing a fourth ’shift’ parameter $\alpha$ to the regular case. This additional degree of freedom allows for the optimization of maneuvers, significantly reducing energy consumption and providing a mechanism to manipulate the gimbal lock condition. While similar optimization strategies have been explored for two-axis sequences of arbitrary length [16,17], this work focuses on an explicit three-axis, four-factor approach, which provides a closed-form expression for the cost function and certainly contains the two-axes setting as a particular case. By utilizing this framework, we can relax the traditional Davenport condition ${\gamma }_{12}+{\gamma }_{23}\ge \pi$, where ${\gamma }_{ij}$ denotes the minimal angle between the rotation axes, to a broader geometric constraint ${\gamma }_{12}+{\gamma }_{23}+{\gamma }_{31}\ge \pi$ allowing for the design of attitude control mechanisms with non-orthogonal gimbals [23,24] that might turn useful in robotics and aircraft engineering.

The text is organized as follows: after a brief preliminary section, the construction for pure rotations is discussed in detail, as derived in [18]. Then, the extension to generic rigid motions is carried out in Section 4, where some issues like the order of finite generations and classical results in spherical geometry, prolonged to the dual sphere, are also discussed. Numerical examples and figures are also provided wherever considered necessary.

2. Generalized Euler Decompositions with Projective Quaternions

The famous Rodrigues’ vector is best understood as a projective quaternion ${\displaystyle \mathbf{c}=\frac{\zeta }{{\zeta }_0}}$, where we can normalize the homogeneous coordinates for convenience

$$\zeta =({\zeta }_0,\zeta )={\zeta }_0+{\zeta }_1\mathbf{i}+{\zeta }_2\mathbf{j}+{\zeta }_3\mathbf{k}\in \mathbb{H},{\parallel \zeta \parallel }^2={\zeta }_0^2+{\zeta }^2=1$$

so that $\zeta$ parametrizes the spin group $\mathrm{SU}\left(2\right)\cong {\mathbb{S}}^3$. This yields a natural representation $\mathrm{SO}\left(3\right)\cong {\mathbb{RP}}^3$ devoid of singularities. Note that one may also invoke Euler’s trigonometric substitution and express $\mathbf{c}=\tau \mathbf{n}$ where ${\displaystyle \tau =tan\frac{\varphi }{2}\in {\mathbb{RP}}^1}$ is usually referred to as the scalar parameter ($\varphi$ being the angle of rotation) and $\mathbf{n}\in {\mathbb{S}}^2$ is the unit normal to the rotation plane.

That way we can express the rotation matrix entries in terms of rational functions only

$$\mathcal{R}\left(\mathbf{c}\right)={\displaystyle \frac{1-{\mathbf{c}}^2+2\mathbf{c}{\mathbf{c}}^t+2{\mathbf{c}}^{\times }}{1+{\mathbf{c}}^2}}={\displaystyle \frac{1+{\mathbf{c}}^{\times }}{1-{\mathbf{c}}^{\times }}}=\mathrm{Cay}\left({\mathbf{c}}^{\times }\right)$$

(1)

vie the well-known Cayley map. Here $\mathbf{c}{\mathbf{c}}^t$ stands for the tensor (dyadic) product and ${\mathbf{c}}^{\times }$ denotes the adjoined, i.e., ${\mathbf{a}}^{\times }\mathbf{b}=\mathbf{a}\times \mathbf{b}$. Another advantage is the simple composition law

$$\langle {\mathbf{c}}_2,{\mathbf{c}}_1\rangle ={\displaystyle \frac{{\mathbf{c}}_2+{\mathbf{c}}_1+{\mathbf{c}}_2\times {\mathbf{c}}_1}{1-{\mathbf{c}}_2·{\mathbf{c}}_1}}\iff \mathcal{R}\left({\mathbf{c}}_2\right)\mathcal{R}\left({\mathbf{c}}_1\right)=\mathcal{R}\left(\langle {\mathbf{c}}_2,{\mathbf{c}}_1\rangle \right)$$

(2)

(here $\mathbf{a}·\mathbf{b}$ stands for the dot product) that is just the projective version of quaternion multiplication. It appears to be far more efficient compared to the usual matrix formalism.

The above description allows us to easily derive (see [4] for details) explicit solutions to the generalized Euler decomposition problem $\mathcal{R}\left(\mathbf{c}\right)=\mathcal{R}\left({\mathbf{c}}_3\right)\mathcal{R}\left({\mathbf{c}}_2\right)\mathcal{R}\left({\mathbf{c}}_1\right)$ with initially given rotation axes ${\mathbf{a}}_i\in {\mathbb{S}}^2$, such that ${\mathbf{c}}_i={\tau }_i{\mathbf{a}}_i$ where ${\displaystyle {\tau }_i=tan\frac{{\varphi }_i}{2}}$. More precisely, we have

$${\tau }_i^{\pm }={\sigma }_i{\displaystyle \frac{{\omega }_i\mp \sqrt{\Delta }}{{\omega }_i^2-\Delta }}={\displaystyle \frac{{\sigma }_i}{{\omega }_i\pm \sqrt{\Delta }}},{\sigma }_i={\epsilon }_{ijk}(g_{jk}-r_{jk}),j>k$$

(3)

where ${\epsilon }_{ijk}$ denotes the Levi-Civita symbol. Furthermore, we make use of the notation

$$g_{ij}={\mathbf{a}}_i·{\mathbf{a}}_j,r_{ij}={\widehat{\mathbf{c}}}_i·\mathcal{R}\left(\mathbf{c}\right){\mathbf{a}}_j$$

as well as

$${\omega }_1={\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathcal{R}}^t\left(\mathbf{c}\right){\mathbf{a}}_3,{\omega }_2={\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathbf{a}}_3,{\omega }_3=\mathcal{R}\left(\mathbf{c}\right){\mathbf{a}}_1·{\mathbf{a}}_2\times {\mathbf{a}}_3$$

and the necessary and sufficient condition for the existence of real solutions

$$\Delta =\left|\begin{array}{ccc}1& g_{12}& r_{31}\\ g_{21}& 1& g_{23}\\ r_{31}& g_{32}& 1\end{array}\right|\ge 0.$$

(4)

Whenever any of the conditions $r_{ji}=g_{ji}$ ($i\ne j$) holds, one may decompose into a pair of rotations about the i-th and the j-th axis as $\mathcal{R}\left(\mathbf{c}\right)=\mathcal{R}\left({\tau }_j{\mathbf{a}}_j\right)\mathcal{R}\left({\tau }_i{\mathbf{a}}_i\right)$, e.g., for $i=1,j=2$:

$${\tau }_1={\displaystyle \frac{r_{22}-1}{{\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathcal{R}}^t\left(\mathbf{c}\right){\mathbf{a}}_2}},{\tau }_2={\displaystyle \frac{r_{11}-1}{\mathcal{R}\left(\mathbf{c}\right){\mathbf{a}}_1·{\mathbf{a}}_1\times {\mathbf{a}}_2}}·$$

(5)

These expressions apply also to the singular gimbal lock setting

$${\mathbf{a}}_3=\pm \mathcal{R}{\mathbf{a}}_1\iff r_{31}=\pm 1$$

(6)

in which the first equality in (5) yields the scalar parameter of ${\varphi }_1\pm {\varphi }_3$ (see also [5] and [4]).

3. The Shift Parameter Construction

This section explains the results obtained in [18] quite thoroughly, since the derivation in the dual quaternion setting follows the same pattern. Let us consider decompositions into four consecutive rotations about three given axes—for the one, which takes part twice, we have an additional ’shift’ angle $\alpha$, that can be varied for the sake of optimization. Denoting the rotation by an angle $\phi$ about ${\mathbf{a}}_i$ with ${\mathcal{R}}_1\left(\phi \right)$, we can write for example

$$\mathcal{R}={\mathcal{R}}_3\left(\psi \right){\mathcal{R}}_1\left(\alpha \right){\mathcal{R}}_2\left(\vartheta \right){\mathcal{R}}_1\left(\phi \right)={\mathcal{R}}_3\left(\psi \right){\mathcal{R}}_2^{\prime }\left(\vartheta \right){\mathcal{R}}_1(\phi +\alpha )$$

(7)

where ${\mathcal{R}}_2^{\prime }$ stands for a rotation about the axis ${\mathbf{a}}_2^{\prime }={\mathcal{R}}_1\left(\alpha \right){\mathbf{a}}_2$ (see for example [5]). The next step is to use the additional fourth parameter $\alpha$ to guarantee the necessary and sufficient condition (4) for the corresponding decomposition, which in this setting takes the form

$${\Delta }^{\prime }=\left|\begin{array}{ccc}1& g_{12}& r_{31}\\ g_{12}& 1& g_{23}^{\prime }\\ r_{31}& g_{23}^{\prime }& 1\end{array}\right|\ge 0$$

(8)

where one may easily express the $\alpha$-dependent component as

$$g_{23}^{\prime }\left(\alpha \right)={\mathbf{a}}_2^{\prime }·{\mathbf{a}}_3=g_{12}{g}_{13}+g_{1[1}{g}_{2]3}cos\alpha +{\omega }_{2}sin\alpha a_{[i}{b}_{j]}=a_i{b}_j-a_j{b}_i.$$

On the other hand, (8) may be written in the equivalent form

$$g_{12}{r}_{31}-\sqrt{(1-g_{12}^2)(1-r_{31}^2)}\le g_{23}^{\prime }\le g_{12}{r}_{31}+\sqrt{(1-g_{12}^2)(1-r_{31}^2)}.$$

Using basic trigonometry one may simplify the expression for $g_{23}^{\prime }$ as

$$g_{23}^{\prime }\left(\alpha \right)=Acos(\alpha -{\alpha }_0)+g_{12}{g}_{13}$$

with

$$A=\sqrt{{\left(g_{1[1}{g}_{2]3}\right)}^2+{\omega }_2^2},{\alpha }_0=arctan\left({\displaystyle \frac{{\omega }_2}{g_{1[1}{g}_{2]3}}}\right)·$$

This allows for writing the above inequalities in the simpler form

$$cos({\gamma }_{12}+{\tilde{\gamma }}_{31})-g_{12}{g}_{13}\le Acos(\alpha -{\alpha }_0)\le cos({\gamma }_{12}-{\tilde{\gamma }}_{31})-g_{12}{g}_{13}$$

where we denote, assuming for the relative angles ${\displaystyle {\gamma }_{ij}\in \left[0,\frac{\pi }{2}\right]}$, while ${\tilde{\gamma }}_{ij}\in [0,\pi ]$

$${\gamma }_{ij}=arccos{g}_{ij}=\measuredangle ({\mathbf{a}}_i,{\mathbf{a}}_j),{\tilde{\gamma }}_{ij}=arccos{r}_{ij}=\measuredangle ({\widehat{\mathbf{c}}}_i,\mathcal{R}{\mathbf{a}}_j).$$

Since we have the restriction $|cos(\alpha -{\alpha }_0)|\le 1$, real solutions are guaranteed if and only if

$$cos({\gamma }_{12}-{\tilde{\gamma }}_{31})-g_{12}{g}_{13}\ge -A,cos({\gamma }_{12}+{\tilde{\gamma }}_{31})-g_{12}{g}_{13}\le A.$$

On the other hand, one has

$${\displaystyle {\omega }_2^2=det\left\{g_{ij}\right\}\Rightarrow A=\sqrt{(1-g_{12}^2)(1-g_{13}^2)}}$$

so the above conditions may be written also as

$$cos({\gamma }_{12}-{\tilde{\gamma }}_{31})\ge cos({\gamma }_{12}+{\gamma }_{31}),cos({\gamma }_{12}+{\tilde{\gamma }}_{31})\le cos({\gamma }_{12}-{\gamma }_{31})$$

(9)

which finally yields

$$|{\gamma }_{12}-{\gamma }_{13}|-{\gamma }_{12}\le {\tilde{\gamma }}_{31}\le 2{\gamma }_{12}+{\gamma }_{13}.$$

(10)

The interval determined in this way is typically larger compared to

$$|{\gamma }_{12}-{\gamma }_{23}|\le {\tilde{\gamma }}_{31}\le {\gamma }_{12}+{\gamma }_{23}$$

(11)

that guarantees decomposition into three consecutive rotations about the ${\mathbf{a}}_i$’s (see [4]). Moreover, since the extremal values of ${\tilde{\gamma }}_{31}$ correspond to $\mathbf{n}\perp {\mathbf{a}}_{1,3}$, one obviously has

$${\gamma }_{31}-\left|\varphi \right|\le {\tilde{\gamma }}_{31}\le {\gamma }_{31}+\left|\varphi \right|$$

(12)

and the sufficient condition for the compound rotation’s angle is

$$\left|\varphi \right|\le 2{\gamma }_{12}.$$

(13)

Note that as long as (10) holds, one may choose an arbitrary

$$\alpha \in \left[{\alpha }_0+arccos{\chi }^{-},{\alpha }_0+arccos{\chi }^{+}\right]$$

(14)

with the notation

$${\displaystyle {\chi }^{\pm }=A^{-1}\left[cos({\gamma }_{12}\pm {\tilde{\gamma }}_{31})-g_{12}{g}_{13}\right].}$$

The endpoints of the interval (14), for example, yield ${\Delta }^{\prime }=0$ and hence, rational expressions for the ${\tau }_k$’s. More generally, denoting (for other values of the indices $g_{ij}^{\prime }=g_{ij}$ and $r_{ij}^{\prime }=r_{ij}$)

$$g_{23}^{\prime }={\mathbf{a}}_2^{\prime }·{\mathbf{a}}_3,r_{21}^{\prime }={\mathbf{a}}_2^{\prime }·\mathcal{R}{\mathbf{a}}_1,r_{32}^{\prime }={\mathbf{a}}_3·\mathcal{R}{\mathbf{a}}_2^{\prime }$$

as well as

$${\omega }_1^{\prime }={\mathbf{a}}_1\times {\mathbf{a}}_2^{\prime }·{\mathcal{R}}^t{\mathbf{a}}_3,{\omega }_2^{\prime }={\mathbf{a}}_1\times {\mathbf{a}}_2^{\prime }·{\mathbf{a}}_3,{\omega }_3^{\prime }=\mathcal{R}{\mathbf{a}}_1·{\mathbf{a}}_2^{\prime }\times {\mathbf{a}}_3$$

one obtains the explicit solutions (3) in the form

$${\tau }_i^{\pm }={\displaystyle \frac{{\sigma }_i^{\prime }}{{\omega }_i^{\prime }\pm \sqrt{{\Delta }^{\prime }}}},{\sigma }_i^{\prime }={\epsilon }_{ijk}(g_{jk}^{\prime }-r_{jk}^{\prime }),j>k$$

(15)

which provides the remaining rotation angles in the decomposition

$${\phi }^{\pm }=2arctan{\tau }_1^{\pm }-\alpha ,{\vartheta }^{\pm }=2arctan{\tau }_2^{\pm },{\psi }^{\pm }=2arctan{\tau }_3^{\pm }.$$

Following a similar approach, one may also consider the decomposition

$$\mathcal{R}={\mathcal{R}}_2\left(\alpha \right){\mathcal{R}}_3\left(\psi \right){\mathcal{R}}_2\left(\vartheta \right){\mathcal{R}}_1\left(\phi \right)={\mathcal{R}}_3^{\prime }\left(\psi \right){\mathcal{R}}_2(\vartheta +\alpha ){\mathcal{R}}_1\left(\phi \right)$$

(16)

with ${\mathbf{a}}_3^{\prime }={\mathcal{R}}_2\left(\alpha \right){\mathbf{a}}_3$ and the corresponding discriminant condition

$${\Delta }^{\prime }=\left|\begin{array}{ccc}1& g_{12}& r_{31}^{\prime }\\ g_{12}& 1& g_{23}\\ r_{31}^{\prime }& g_{23}& 1\end{array}\right|\ge 0$$

(17)

where the $\alpha$-dependent component $r_{31}^{\prime }$ may be written explicitly in the form

$$r_{31}^{\prime }\left(\alpha \right)={\mathbf{a}}_3^{\prime }·\mathcal{R}{\mathbf{a}}_1=g_{2[2}{r}_{3]1}cos\alpha +{\omega }_{3}sin\alpha +g_{23}{r}_{21}$$

and further simplified as

$$r_{31}^{\prime }\left(\alpha \right)=\tilde{A}cos(\alpha -{\alpha }_0)+g_{23}{r}_{21}$$

$$\tilde{A}=\sqrt{{\left(g_{2[2}{r}_{3]1}\right)}^2+{\omega }_3^2},{\alpha }_0=arctan\left({\displaystyle \frac{{\omega }_3}{g_{2[2}{r}_{3]1}}}\right)·$$

This yields the inequalities

$$cos({\gamma }_{12}+{\gamma }_{23})-g_{23}{r}_{21}\le \tilde{A}cos(\alpha -{\alpha }_0)\le cos({\gamma }_{12}-{\gamma }_{23})-g_{23}{r}_{21}$$

which allow real solutions for $\alpha$ as long as

$$cos({\gamma }_{12}+{\gamma }_{23})-g_{23}{r}_{21}\le \tilde{A},cos({\gamma }_{12}-{\gamma }_{23})-g_{23}{r}_{21}\ge -\tilde{A}.$$

Furthermore, taking into account that

$${\omega }_3^2=\left|\begin{array}{ccc}1& g_{23}& r_{21}\\ g_{23}& 1& r_{31}\\ r_{21}& r_{31}& 1\end{array}\right|\Rightarrow \tilde{A}=\sqrt{(1-g_{23}^2)(1-r_{21}^2)}$$

one ends up with

$$cos({\gamma }_{23}-{\tilde{\gamma }}_{21})\ge cos({\gamma }_{12}+{\gamma }_{23}),cos({\gamma }_{23}+{\tilde{\gamma }}_{21})\le cos({\gamma }_{12}-{\gamma }_{23})$$

(18)

which yields for ${\tilde{\gamma }}_{21}$ the interval

$$|{\gamma }_{23}-{\gamma }_{12}|-{\gamma }_{23}\le {\tilde{\gamma }}_{21}\le 2{\gamma }_{23}+{\gamma }_{12}$$

(19)

and with the aid of formula (12), one obtains a sufficient condition for $\varphi$ in the form

$$\left|\varphi \right|\le 2{\gamma }_{23}.$$

(20)

In the symmetric case ${\gamma }_{12}={\gamma }_{23}=\gamma$ one has ${\tilde{\gamma }}_{21}\le 3\gamma$, while for $r_{21}=\pm 1$, (14) yields either no solution or no restriction for $\alpha$, the condition for the latter being ${\gamma }_{12}\le 2{\gamma }_{23}$ and ${\gamma }_{12}\ge \pi -2{\gamma }_{23}$, respectively. If (19) holds and $\tilde{A}\ne 0$, then $\alpha$ is arbitrary within (14), where

$${\displaystyle {\chi }^{\pm }={\tilde{A}}^{-1}\left[cos({\gamma }_{12}\pm {\gamma }_{23})-g_{23}{r}_{31}\right]}$$

and calculate (for all other values of i and j we have $g_{ij}^{\prime }=g_{ij}$, respectively $r_{ij}^{\prime }=r_{ij}$)

$$g_{31}^{\prime }={\mathbf{a}}_3^{\prime }·{\mathbf{a}}_1,r_{32}^{\prime }={\mathbf{a}}_3^{\prime }·\mathcal{R}{\mathbf{a}}_2,r_{31}^{\prime }={\mathbf{a}}_3^{\prime }·\mathcal{R}{\mathbf{a}}_1$$

$${\omega }_1^{\prime }={\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathcal{R}}^t{\mathbf{a}}_3^{\prime },{\omega }_2^{\prime }={\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathbf{a}}_3^{\prime },{\omega }_3^{\prime }=\mathcal{R}{\mathbf{a}}_1·{\mathbf{a}}_2\times {\mathbf{a}}_3^{\prime }$$

in order to obtain the solutions based on (3) in the form

$${\tau }_i^{\pm }={\displaystyle \frac{{\sigma }_i^{\prime }}{{\omega }_i^{\prime }\pm \sqrt{{\Delta }^{\prime }}}},{\sigma }_i^{\prime }={\epsilon }_{ijk}(g_{jk}^{\prime }-r_{jk}^{\prime }),j>k$$

(21)

which finally provides the remaining three rotation angles

$${\phi }^{\pm }=2arctan{\tau }_1^{\pm },{\vartheta }^{\pm }=2arctan{\tau }_2^{\pm }-\alpha ,{\psi }^{\pm }=2arctan{\tau }_3^{\pm }.$$

Note also that (16) is, in some sense, dual to (7), since former is equivalent to

$${\mathcal{R}}^t={\mathcal{R}}_1(-\phi ){\mathcal{R}}_3^{\prime \prime }(-\psi ){\mathcal{R}}_2(-\vartheta -\alpha )$$

(22)

where ${\mathcal{R}}_3^{\prime \prime }$ is a rotation about ${\mathbf{a}}_3^{\prime \prime }={\mathcal{R}}_2(-\vartheta ){\mathbf{a}}_3$, so the shift angle in this case is $\vartheta$ (see [18]):

$$\vartheta \in \left[{\vartheta }_0+arccos{\chi }^{-},{\vartheta }_0+arccos{\chi }^{+}\right]$$

(23)

$${\chi }^{\pm }={\displaystyle \frac{cos({\gamma }_{23}\pm {\tilde{\gamma }}_{21})-g_{12}{g}_{23}}{sin{\gamma }_{12}sin{\gamma }_{23}}},{\vartheta }_0=arctan\left({\displaystyle \frac{{\omega }_2}{g_{1[2}{g}_{3]2}}}\right)·$$

It is straightforward to see that the above expression takes real values as long as

$$cos({\gamma }_{23}+{\tilde{\gamma }}_{21})-g_{12}{g}_{23}\le \tilde{A}cos(\vartheta -{\vartheta }_0)\le cos({\gamma }_{23}-{\tilde{\gamma }}_{21})-g_{12}{g}_{23}$$

$${\displaystyle \tilde{A}=\sqrt{{\left(g_{1[2}{g}_{3]2}\right)}^2+{\omega }_2^2}=sin{\gamma }_{12}sin{\gamma }_{23},{\vartheta }_0=arctan\left(\frac{{\omega }_2}{g_{1[2}{g}_{3]2}}\right)}$$

which yields (19), and if $\tilde{A}\ne 0$, one may choose arbitrary $\vartheta$ in the interval (23) and calculate

$$g_{21}^{\prime \prime }=g_{23},g_{31}^{\prime \prime }=g_{12},g_{32}^{\prime \prime }={\mathbf{a}}_3^{\prime \prime }·{\mathbf{a}}_1,r_{31}^{\prime \prime }=r_{21},r_{32}^{\prime \prime }={\mathbf{a}}_3^{\prime \prime }·\mathcal{R}{\mathbf{a}}_1,r_{21}^{\prime \prime }={\mathbf{a}}_2·\mathcal{R}{\mathbf{a}}_3^{\prime \prime }$$

$${\omega }_1^{\prime \prime }={\mathbf{a}}_2\times {\mathbf{a}}_3^{\prime \prime }·\mathcal{R}{\mathbf{a}}_1,{\omega }_2^{\prime \prime }={\mathbf{a}}_1\times {\mathbf{a}}_2·{\mathbf{a}}_3^{\prime \prime },{\omega }_3^{\prime \prime }={\mathcal{R}}^t{\mathbf{a}}_2·{\mathbf{a}}_3^{\prime \prime }\times {\mathbf{a}}_1.$$

Finally, ${\Delta }^{\prime \prime }$ is given by (4) with $g_{ij}^{\prime \prime }$ and $r_{31}^{\prime \prime }$ instead of $g_{ij}$ and $r_{31}$, respectively, and one has

$${\tau }_i^{\pm }={\displaystyle \frac{{\sigma }_i^{\prime \prime }}{{\omega }_i^{\prime \prime }\pm \sqrt{{\Delta }^{\prime \prime }}}},{\sigma }_i^{\prime \prime }={\epsilon }_{ijk}(g_{jk}^{\prime \prime }-r_{jk}^{\prime \prime }),j>k$$

(24)

from formula (3), which yields the remaining three rotation angles

$${\phi }^{\pm }=-2arctan{\tau }_3^{\pm },{\psi }^{\pm }=-2arctan{\tau }_2^{\pm },{\alpha }^{\pm }=-\vartheta -2arctan{\tau }_1^{\pm }.$$

Next, let us consider decompositions of the type

$$\mathcal{R}=\mathcal{R}\left(\mathbf{c}\right)={\mathcal{R}}_1\left(\alpha \right){\mathcal{R}}_3\left(\psi \right){\mathcal{R}}_2\left(\vartheta \right){\mathcal{R}}_1\left(\phi \right)$$

(25)

in which we conjugate the compound rotation rather than a particular factor, namely

$$\tilde{\mathcal{R}}=\mathcal{R}\left({\mathcal{R}}_1^t\left(\alpha \right)\mathbf{c}\right)={\mathcal{R}}_3\left(\psi \right){\mathcal{R}}_2\left(\vartheta \right){\mathcal{R}}_1(\phi +\alpha ).$$

Like in the previous cases, we express the $\alpha$-dependent parameter

$$r_{31}^{\prime }\left(\alpha \right)={\mathcal{R}}_1\left(\alpha \right){\mathbf{a}}_3·\mathcal{R}{\mathbf{a}}_1=g_{1[1}{r}_{3]1}cos\alpha +{\tilde{\omega }}_{2}sin\alpha +g_{13}{r}_{11}$$

using the notation ${\tilde{\omega }}_2={\mathbf{a}}_1\times {\mathbf{a}}_3·\mathcal{R}{\mathbf{a}}_1$, in the form

$$r_{31}^{\prime }\left(\alpha \right)=\tilde{A}cos(\alpha -{\alpha }_0)+g_{13}{r}_{11}$$

$$\tilde{A}=\sqrt{{\left(g_{1[1}{r}_{3]1}\right)}^2+{\tilde{\omega }}_2^2},{\alpha }_0=arctan\left({\displaystyle \frac{{\tilde{\omega }}_2}{g_{1[1}{r}_{3]1}}}\right)·$$

Then, the necessary and sufficient condition

$${\Delta }^{\prime }=\left|\begin{array}{ccc}1& g_{12}& r_{31}^{\prime }\\ g_{12}& 1& g_{23}\\ r_{31}^{\prime }& g_{23}& 1\end{array}\right|\ge 0$$

(26)

may be expressed in terms of the shift angle $\alpha$ as

$$cos({\gamma }_{12}+{\gamma }_{23})-g_{13}{r}_{11}\le \tilde{A}cos(\alpha -{\alpha }_0)\le cos({\gamma }_{12}-{\gamma }_{23})-g_{13}{r}_{11}.$$

On the other hand, one has

$${\tilde{\omega }}_2^2=\left|\begin{array}{ccc}1& g_{13}& r_{11}\\ g_{13}& 1& r_{31}\\ r_{11}& r_{31}& 1\end{array}\right|\Rightarrow \tilde{A}=\sqrt{(1-g_{13}^2)(1-r_{11}^2)}$$

so real solutions for $\alpha$ exist as long as

$$cos({\gamma }_{13}-{\tilde{\gamma }}_{11})\ge cos({\gamma }_{12}+{\gamma }_{23}),cos({\gamma }_{13}+{\tilde{\gamma }}_{11})\le cos({\gamma }_{12}-{\gamma }_{23})$$

(27)

that yields for ${\tilde{\gamma }}_{11}$ the interval

$$|{\gamma }_{23}-{\gamma }_{12}|-{\gamma }_{13}\le {\tilde{\gamma }}_{11}\le {\gamma }_{12}+{\gamma }_{23}+{\gamma }_{13}.$$

(28)

This time the lower bound is trivial since ${\tilde{\gamma }}_{11}$ is non-negative, while the expression on the left is non-positive due to the triangle inequality. Thus the sufficient condition for the angle

$$\left|\varphi \right|\le {\gamma }_{12}+{\gamma }_{23}+{\gamma }_{13}.$$

(29)

In the generic case we choose a value of $\alpha$ in the non-empty interval (14), where

$${\displaystyle {\chi }^{\pm }={\tilde{A}}^{-1}\left[cos({\gamma }_{12}\pm {\gamma }_{23})-g_{13}{r}_{11}\right]}$$

assuming $\tilde{A}\ne 0$ ($\tilde{A}=0$ does not impose restriction on $\alpha$) and use it to obtain the shifted rotation $\tilde{\mathcal{R}}=\mathcal{R}\left({\mathcal{R}}_1^t\left(\alpha \right)\mathbf{c}\right)={\mathcal{R}}_1^t\left(\alpha \right)\mathcal{R}\left(\mathbf{c}\right){\mathcal{R}}_1\left(\alpha \right)$ that is to replace $\mathcal{R}\left(\mathbf{c}\right)$ in (3). Then, calculate

$$r_{ij}^{\prime }={\widehat{\mathbf{c}}}_i^{\prime }·\mathcal{R}\left(\mathbf{c}\right){\mathbf{a}}_j^{\prime },{\widehat{\mathbf{c}}}_i^{\prime }={\mathcal{R}}_1\left(\alpha \right){\mathbf{a}}_i$$

with ${\mathbf{a}}_i^{\prime }={\mathcal{R}}_1\left(\alpha \right){\mathbf{a}}_i$, as well as

$${\omega }_1^{\prime }={\mathbf{a}}_1\times {\mathbf{a}}_2^{\prime }·{\mathcal{R}}^t\left(\mathbf{c}\right){\mathbf{a}}_3^{\prime },{\omega }_2^{\prime }={\omega }_2,{\omega }_3^{\prime }=\mathcal{R}\left(\mathbf{c}\right){\mathbf{a}}_1·{\mathbf{a}}_2^{\prime }\times {\mathbf{a}}_3^{\prime }$$

and thus express

$${\tau }_i^{\pm }={\displaystyle \frac{{\sigma }_i^{\prime }}{{\omega }_i^{\prime }\pm \sqrt{{\Delta }^{\prime }}}},{\sigma }_i^{\prime }={\epsilon }_{ijk}(g_{jk}-r_{jk}^{\prime }),j>k$$

(30)

where ${\Delta }^{\prime }$ is given by formula (26), finally arriving at

$${\psi }^{\pm }=2arctan{\tau }_3^{\pm },{\vartheta }^{\pm }=2arctan{\tau }_2^{\pm },{\phi }^{\pm }=2arctan{\tau }_1^{\pm }-\alpha .$$

3.1. A Brief Note on Gimbal Lock Control

Gimbal lock (6) is a zero-measure set singularity so avoiding it is easy, as long as the corresponding interval for the shift parameter is non-degenerate, i.e., contains more than one point. Invoking it, however, is less trivial: consider for example (16), where the condition for the shifted third axis ${\mathbf{a}}_3^{\prime }={\mathcal{R}}_2\left(\alpha \right){\mathbf{a}}_3=\pm \mathcal{R}{\mathbf{a}}_1$ demands that (see also [7])

$${\mathbf{a}}_2\in \mathrm{Span}\{{\mathbf{a}}_3\pm \mathcal{R}{\mathbf{a}}_1,{\mathbf{a}}_3\times \mathcal{R}{\mathbf{a}}_1\}\iff (1\pm r_{31})(r_{21}\mp g_{23})=0.$$

(31)

The non-trivial solutions $r_{21}=\pm g_{23}$ impose the conditions ${\gamma }_{12}={\gamma }_{23}$ and respectively ${\gamma }_{12}+{\gamma }_{23}=\pi$ (see [18]), allowing for factorizations into pairs, since $\mathcal{R}={\mathcal{R}}_2{\mathcal{R}}_1\iff r_{21}=g_{21}$. where ${\mathcal{R}}_2^{\prime }$ stands for a rotation about the unit vector ${\mathbf{a}}_2^{\prime }={\mathcal{R}}_1(-\phi ){\mathbf{a}}_2$, so $\phi$ plays the role of a shift angle. Resolving the factorization like before, we end up with a gimbal lock condition in the form $(1\pm r_{32})(r_{31}\mp g_{12})=0$ with non-trivial solutions $r_{31}=\pm g_{21}$ demanding ${\gamma }_{12}={\gamma }_{13}$ or ${\gamma }_{12}+{\gamma }_{13}=\pi$, respectively. In both cases, however, we also have $r_{31}=g_{31}$ which guarantees the factorization $\mathcal{R}={\mathcal{R}}_3{\mathcal{R}}_1$.

Finally, the gimbal lock condition in (25) may be written as $(1\pm r_{31})(r_{11}\mp g_{31})=0$ and obviously has a non-trivial solution in the form $r_{11}=\pm g_{31}$. Since $r_{31}^{\prime }=\pm 1$ in this case is equivalent to $cos(\alpha -{\alpha }_0)=\pm 1$, it demands ${\gamma }_{12}={\gamma }_{23}$ or ${\gamma }_{12}+{\gamma }_{23}=\pi$, respectively. In [18] we conclude that avoiding gimbal lock this way is always possible as long as there are two distinct angles ${\gamma }_{ij}$, whose sum is non less than ${\displaystyle \frac{\pi }{2}}$, while inflicting it may be done only in one of the two cases: (i.) $r_{ij}=g_{ij}=\pm g_{jk}=0$ that yields a decomposition into two factors, and (ii.) $r_{ii}=\pm g_{ij}$ with $g_{ik}=\pm g_{jk}$, where the values of $i,j$ and k are assumed to be different (allowing for permutations of axes). Here we skip some calculations, for brevity.

3.2. Optimization

Here we show how the additional parameter $\alpha$ is used to minimize the cost function

$${\mathcal{E}}_{\gamma }=\left|\phi \right|+\left|\vartheta \right|+\left|\psi \right|+\left|\alpha \right|$$

(32)

with a few examples. Consider (7) and let $\left\{{\mathbf{a}}_i\right\}$ be the standard basis, i.e., $g_{ij}={\delta }_{ij}$ and $r_{ij}$ are simply the matrix entries of $\mathcal{R}$. Since the $\alpha -$rotation of ${\mathbf{a}}_2$ takes place in the $YOZ$ plane, it can be written as ${\mathbf{a}}_2^{\prime }=cos\alpha {\mathbf{a}}_2+sin\alpha {\mathbf{a}}_3$ that results in the expressions

$$g_{32}^{\prime }=sin\alpha ,r_{21}^{\prime }=r_{21}cos\alpha +r_{31}sin\alpha ,r_{32}^{\prime }=r_{32}cos\alpha +r_{33}sin\alpha$$

$${\omega }_1^{\prime }=r_{33}cos\alpha -r_{32}sin\alpha ,{\omega }_2^{\prime }=cos\alpha ,{\omega }_3^{\prime }=r_{11}cos\alpha$$

and ${\displaystyle {\Delta }^{\prime }={cos}^2\alpha -r_{31}^2\Rightarrow \left|\alpha \right|\le arcsin\sqrt{1-r_{31}^2}}$. Substituting in formula (15) then yields

$$\begin{array}{c}\hfill {\tau }_1^{\pm }={\displaystyle \frac{r_{32}+(r_{33}-1)tan\alpha }{{\displaystyle r_{33}-r_{32}tan\alpha \pm \sqrt{1-r_{31}^2{sec}^2\alpha }}}}=tan{\displaystyle \frac{\alpha +\phi }{2}}\\ & \\ \hfill {\tau }_2^{\pm }=-{\displaystyle \frac{r_{31}sec\alpha }{{\displaystyle 1\pm \sqrt{1-r_{31}^2{sec}^2\alpha }}}}=tan{\displaystyle \frac{\vartheta }{2}},{\tau }_3^{\pm }={\displaystyle \frac{r_{21}+r_{31}tan\alpha }{{\displaystyle r_{11}\pm \sqrt{1-r_{31}^2{sec}^2\alpha }}}}=tan{\displaystyle \frac{\psi }{2}}\end{array}$$

(33)

and the cost function is finally expressed in the form

$${\displaystyle {\mathcal{E}}_{\gamma }\left(\alpha \right)=\left|\alpha \right|+|2arctan{\tau }_1-\alpha |+|2arctan{\tau }_2|+|2arctan{\tau }_3|.}$$

The gimbal lock singularity (6) persists for all values of $\alpha$, but (8) demands $\alpha =0$ or $\alpha =\pi$, so the free parameter is now $\psi \in [-\pi ,\pi ]$ and minimize ${\mathcal{E}}_{\gamma }\left(\psi \right)$ instead. Figure 1a illustrates (32) for a rotation by an angle $\varphi =-{120}^{\circ }$ about the vector $\mathbf{n}\sim {(3,4,5)}^t$. Its minimum ${\mathcal{E}}_{\gamma }\left({\alpha }_0\right)\approx 179.{81}^{\circ }$ is reached at ${\alpha }_0\approx 13.{49}^{\circ }$, that is about $36.49\%$ more efficient compared to the standard Bryan setting ${\mathcal{E}}_{\gamma }\left(0\right)\approx 245.{31}^{\circ }$, and thus the optimal angles are

$$\phi \approx -0.{16}^{\circ },\vartheta \approx -75.{14}^{\circ },\alpha \approx 13.{49}^{\circ },\psi \approx -91.{02}^{\circ }.$$

Next , we consider (16) with ${\mathbf{a}}_1={\mathbf{a}}_3$ oriented along the z-coordinate axis and ${\mathbf{a}}_2$ aligned with $OX$, which yields ${\mathbf{a}}_3^{\prime }=cos\alpha {\widehat{\mathbf{e}}}_z-sin\alpha {\widehat{\mathbf{e}}}_y$ and hence, for the scalar parameters

$$\begin{array}{c}\hfill {\tau }_1^{\pm }={\displaystyle \frac{r_{31}-r_{21}tan\alpha }{r_{32}-r_{22}tan\alpha \pm \sqrt{{\Delta }^{\prime }}}}=tan{\displaystyle \frac{\phi }{2}},{\tau }_2^{\pm }={\displaystyle \frac{1-r_{33}+r_{23}tan\alpha }{-tan\alpha \pm \sqrt{{\Delta }^{\prime }}}}=tan{\displaystyle \frac{\vartheta +\alpha }{2}}\\ \hfill {\tau }_3^{\pm }=-{\displaystyle \frac{r_{13}sec\alpha }{r_{23}+r_{33}tan\alpha \mp \sqrt{{\Delta }^{\prime }}}}=tan{\displaystyle \frac{\psi }{2}}\end{array}$$

(34)

where ${\Delta }^{\prime }={sec}^2\alpha -{(r_{33}-r_{23}tan\alpha )}^2$ and by abuse of notation, we let $r_{ij}$ denote the matrix entries of $\mathcal{R}$ in the standard basis. The shift parameter is now unrestricted $\alpha \in [-\pi ,\pi ]$ since $r_{33}^2+r_{23}^2\le 1$. Consider for example a half-turn about ${\displaystyle \mathbf{n}\sim {\left(5,4,3\right)}^t}$ (Figure 1b) with

$${\displaystyle {\mathcal{E}}_{\gamma }\left(\alpha \right)=\left|\alpha \right|+|2arctan{\tau }_1|+|2arctan{\tau }_2-\alpha |+|2arctan{\tau }_3|.}$$

Here $\alpha =0$ is a point of local maximum for ${\mathcal{E}}_{\gamma }\left(\alpha \right)$ and the global minimum at ${\alpha }_0\approx -105.{37}^{\circ }$ yields a total difference of ${\mathcal{E}}_{\gamma }\left(0\right)-{\mathcal{E}}_{\gamma }\left({\alpha }_0\right)\approx 50.{65}^{\circ }$, i.e., about $16.35\%$ decrease in energy consumption, or approximately $19.54\%$ increase in efficiency compared to the standard Euler decomposition. The two points of discontinuity ${\alpha }_1\approx -143.{13}^{\circ }$ and ${\alpha }_2\approx 36.{87}^{\circ }$ are related to the condition $r_{32}^{\prime }=g_{32}^{\prime }$, i.e., $\phi =0$. The optimal sequence then takes the form

$$\phi \approx -52.{24}^{\circ },\vartheta \approx -50.{77}^{\circ },\psi \approx -50.{77}^{\circ },\alpha \approx -105.{37}^{\circ }.$$

Finally, for the $XYZX$ decomposition setting in formula (25) one has

$$\begin{array}{c}\hfill {\tau }_1={\displaystyle \frac{r_{32}-r_{23}{tan}^2\alpha +(r_{33}-r_{22})tan\alpha }{r_{33}+r_{22}{tan}^2\alpha -(r_{32}+r_{23})tan\alpha \pm sec\alpha \sqrt{{\Delta }^{\prime \prime }}}}=tan{\displaystyle \frac{\phi +\alpha }{2}}\\ & \\ \hfill {\tau }_2={\displaystyle \frac{r_{21}tan\alpha -r_{31}}{sec\alpha \pm \sqrt{{\Delta }^{\prime \prime }}}}=tan{\displaystyle \frac{\vartheta }{2}},{\tau }_3={\displaystyle \frac{r_{21}+r_{31}tan\alpha }{r_{11}sec\alpha \pm \sqrt{{\Delta }^{\prime \prime }}}}=tan{\displaystyle \frac{\psi }{2}}\end{array}$$

(35)

with ${\Delta }^{\prime \prime }={sec}^2\alpha -{(r_{31}-r_{21}tan\alpha )}^2$. Since the condition $r_{11}=0$ holds, the endpoints of the interval for $\alpha$ correspond to gimbal lock singularity. Consider again the former rotation discussed above in this section (Figure 1c). The minimum at ${\alpha }_0\approx -74.{97}^{\circ }$ yields about $35.82\%$ increase in efficiency, and the optimal sequence of rotation angles in this case is

$$\phi \approx =2.{39}^{\circ },\vartheta \approx -0.{18}^{\circ },\psi \approx -103.3^{\circ },\alpha \approx -74.8^{\circ }.$$

The cost functions ${\mathcal{E}}_{\gamma }\left(\alpha \right)$ in all examples above imply either spherical or cylindrical symmetry, but in the case of significant anisotropy, this is no longer the case, due to different moments of inertia or external forces. This can be modeled by introducing weights in (32):

$${\mathcal{E}}_{\gamma }={\kappa }_1\left(\right|\phi |+\left|\alpha \right|)+{\kappa }_2\left|\vartheta \right|+{\kappa }_3\left|\psi \right|,{\kappa }_i\in [0,1].$$

Figure 2 illustrates this idea. Note that two of the graphs only udergo a slight deformation, while in Figure 2b we see a phase transition: the maximum at $\alpha =0$ becomes a minimum.

A large part of this section revises an old result by the author, but having the formulas written explicitly here instead of referring to different pages elsewhere should be helpful for a smooth transition to the case of screw motions, considered in the next section. Most of the results follow from the transfer principle, but there are still some subtleties to point out.

4. Dual Quaternions, Screws, and the Transfer Principle

This section is dedicated to implementing the dual extension and applying the transfer principle, which would generalize the above results concerning pure rotations, to arbitrary rigid transformations, interpreted as screw motions in view of the famous Mozzi-Chasles theorem. We begin with the construction of dual numbers $\mathbb{R}\left[\epsilon \right]$ as a central extension to $\mathbb{R}$ with a nilpotent element $\epsilon$ (with ${\epsilon }^2=0$) such that for each $\underline{x}=x+y\epsilon \in \mathbb{R}\left[\epsilon \right]$ we refer to $x=Re\left(\underline{x}\right)$ as the real part of $\underline{x}$ and $y=Du\left(\underline{x}\right)$ as its dual part. We may also define dual conjugate as ${\underline{x}}^{*}=x-y\epsilon$ in analogy with complex numbers, so that $2Re\left(\underline{x}\right)=\underline{x}+{\underline{x}}^{*}$, but retrieving the dual part is less trivial since ${\epsilon }^{-1}$ does not exist. Next, we study the extension to mappings. In particular, for a differentiable function $f\left(x\right)$ we have a finite Taylor series

$$f\left(\underline{x}\right)=f\left(x\right)+f^{\prime }\left(x\right)y\epsilon$$

(36)

due to the defining property of $\epsilon$. This yields for example the trigonometric identities

$$sin(\phi +\epsilon d)=sin\phi +\epsilon dcos\phi ,cos(\phi +\epsilon d)=cos\phi -\epsilon dsin\phi .$$

(37)

Once we have the construction for real numbers, it is rather straightforward to define $\mathbb{C}\left[\epsilon \right]$, $\mathbb{H}\left[\epsilon \right]$, dual vectors and matrices, etc. In particular, the description of 3D Euclidean motions (rigid displacement) uses unit dual quaternions, which constitute the unit dual 3-sphere:

$${\mathbb{S}}^3\left[\epsilon \right]=\{\underline{\mathcal{Q}}\in \mathbb{H}\left[\epsilon \right],|\underline{\mathcal{Q}}{|}^2={\underline{\mathcal{Q}}}^{†}\underline{\mathcal{Q}}=1\}$$

(38)

and are, in essence, pairs of unit and pure quaternions (see [11] for more details)

$$\underline{\mathcal{Q}}=\left(1+{\displaystyle \frac{\epsilon }{2}}t\right)q,{\left|q\right|}^2=1,t^{†}=-t,{\epsilon }^2=0$$

(39)

where q and t represent rotations and translations, respectively. Consider the dual vector

$$\underline{\mathbf{n}}=\mathbf{n}+\epsilon \mathbf{m}\in {\mathbb{S}}^2\left[\epsilon \right]\iff {\mathbf{n}}^2=1,\mathbf{m}\perp \mathbf{n}$$

(40)

and the dual angle $\underline{\phi }=\phi +\epsilon d$ linked to screw displacement d and moment $\mathbf{m}$ defined as

$$d=\mathbf{t}·\mathbf{n},\mathbf{m}={\displaystyle \frac{1}{2}}\left(cot{\displaystyle \frac{\phi }{2}}-\mathbf{n}\times \right){\mathcal{P}}_{\mathbf{n}}^{\perp }\mathbf{t}$$

(41)

with ${\mathcal{P}}_{\mathbf{n}}^{\perp }$ denoting the orthogonal projector associated with $\mathbf{n}$. Hence, the grade projections

$${\langle \underline{\mathcal{Q}}\rangle }_0=q_0-{\displaystyle \frac{\epsilon }{2}}d=cos{\displaystyle \frac{\underline{\phi }}{2}},{\langle \underline{\mathcal{Q}}\rangle }_2=\mathbf{q}+{\displaystyle \frac{\epsilon }{2}}\left(q_0-\mathbf{q}\times \right)\mathbf{t}=\underline{\mathbf{n}}sin{\displaystyle \frac{\underline{\phi }}{2}}$$

(42)

for which we use the correspondence ${\langle q\rangle }_2\iff \mathbf{q}$, formula (39), as well as the identities (37). The famous Mozzi-Chasles theorem asserts that each rigid motion can be represented as a screw displacement, i.e., a rotation and translation with a common axis. This so-called ’screw axis’ is denoted with $\underline{\mathbf{n}}$ above and the dual quaternion construction provides its Plücker coordinates, while the dual angle $\underline{\phi }$ encodes both the rotational and translational component of this motion. We can derive this information from the dual quaternion and conversely, substitute it in formula (42), or use some other familiar representation, e.g., $exp\left(\underline{\phi }{\underline{\mathbf{n}}}^{\times }\right)$ for the explicit form of rigid motions. This analogy is guaranteed by the famous Kotelnikov’s principle of transference (or simply, transfer principle) allowing us to extend all known results for classical quaternions to their dual version, as long as the real part is non-vanishing (otherwise one needs to be concerned with zero divisors). Of course, pure translations may also be incorporated, but with a little bit of caution, as we shall see below.

4.1. The Classical Decomposition of Screw Displacements

We begin with generalizing the older result of Euler-type decompositions to the entire Galilean group, using dual quaternions, and more precisely, their projective version. Given three points on the dual sphere ${\underline{\mathbf{a}}}_k\in {\mathbb{S}}^2\left[\epsilon \right]$ interpreted as oriented screw axes, and one in dual projective space $\underline{\mathbf{c}}\in {\mathbb{RP}}^3\left[\epsilon \right]$ we want to find the values of the scalar parameters ${\underline{\tau }}_k\in {\mathbb{RP}}^1\left[\epsilon \right]$ (if they exist), so that the rigid transformation presented by $\underline{\mathbf{c}}$ is factorized as

$$\underline{\mathbf{c}}=\langle {\underline{\mathbf{c}}}_3,{\underline{\mathbf{c}}}_2,{\underline{\mathbf{c}}}_1\rangle ,{\underline{\mathbf{c}}}_k={\underline{\tau }}_k{\underline{\mathbf{a}}}_k$$

(43)

where we use the associativity of the composition law (2), now extended to the non-homogeneous setting via Kotelnikov’s principle of transference. By the same argument, we can also use formula (1), and thus derive the solutions (3). Note that the condition (4) remains valid only for the real part of the discriminant $\underline{\Delta }=\Delta +\psi \epsilon$ since from (36) we have

$$\sqrt{\Delta +\psi \epsilon }=\sqrt{\Delta }+{\displaystyle \frac{\psi \epsilon }{2\sqrt{\Delta }}}$$

(44)

where $\psi =Du\left(\underline{\Delta }\right)$ can be easily calculated from the components, namely

$$\begin{array}{ccc}\hfill {\underline{\mathbf{a}}}_i={\mathbf{a}}_i+\epsilon {\mathbf{b}}_i& \Rightarrow & {\underline{g}}_{ij}={\underline{\mathbf{a}}}_i·{\underline{\mathbf{a}}}_j={\mathbf{a}}_i·{\mathbf{a}}_j+\epsilon ({\mathbf{a}}_i·{\mathbf{b}}_j+{\mathbf{b}}_i·{\mathbf{a}}_j)=g_{ij}+\epsilon h_{ij}\hfill \\ \hfill \underline{\mathcal{R}}=\mathcal{R}+\epsilon \mathcal{S}& \Rightarrow & {\underline{r}}_{ij}={\mathbf{a}}_i·\mathcal{R}{\mathbf{a}}_j+\epsilon \left({\mathbf{a}}_i·\mathcal{R}{\mathbf{b}}_j+{\mathbf{b}}_i·\mathcal{R}{\mathbf{a}}_j+{\mathbf{a}}_i·\mathcal{S}{\mathbf{a}}_j\right)=r_{ij}+\epsilon s_{ij}\hfill \end{array}$$

ultimately leading to

$$\psi =g_{12}{g}_{23}{s}_{31}+(g_{12}{h}_{23}+h_{12}{g}_{23})r_{31}-2(g_{12}{h}_{12}+g_{23}{h}_{23}+r_{31}{s}_{31}).$$

(45)

Similarly, the dual Rodrigues’ vector is given as

$$\underline{\mathbf{c}}=\mathbf{c}+\epsilon \mathbf{d}=tan{\displaystyle \frac{\varphi }{2}}\left(1+{\displaystyle \frac{\epsilon d}{sin\varphi }}\right)\underline{\mathbf{n}},\mathbf{d}={\displaystyle \frac{d\mathbf{n}+sin\varphi \mathbf{m}}{1+cos\varphi }}$$

(46)

which yields ${\underline{\mathbf{c}}}^2={\mathbf{c}}^2(1+2p\epsilon )$ where $p=dcsc\varphi$ denotes the screw pitch, and respectively

$$\mathcal{S}\left(\underline{\mathbf{c}}\right)=dsin\varphi \left(\mathbf{n}{\mathbf{n}}^t-1\right)+(1-cos\varphi )\left(\mathbf{n}{\mathbf{m}}^t+\mathbf{m}{\mathbf{n}}^t\right)+dcos\varphi {\mathbf{n}}^{\times }+sin\varphi {\mathbf{m}}^{\times }.$$

(47)

Note, however, that the pitch p is ill-defined at both $\varphi =0$ and $\varphi =\pi$, thus it is necessary to clarify the construction for the Rodrigues’ vector itself at these points, unlike in the case of pure rotations, where they represent the trivial element and a half-turn, respectively. At $\varphi =0$ and $d\ne 0$ one has a pure translation, for which the Rodrigues’ vector takes the form

$${\underline{\mathbf{c}}}_i={\displaystyle \frac{\epsilon d_i}{2}}{\mathbf{n}}_i\Rightarrow \langle {\underline{\mathbf{c}}}_1,{\underline{\mathbf{c}}}_2\rangle =\langle {\underline{\mathbf{c}}}_2,{\underline{\mathbf{c}}}_1\rangle ={\displaystyle \frac{\epsilon }{2}}\left(d_1{\mathbf{n}}_1+d_2{\mathbf{n}}_2\right)$$

(48)

while in the case of half-turn screws, although (46) is at infinity, (1) is easily expressed as

$$\underline{\mathcal{R}}\left(\underline{\mathbf{c}}\right)=2\underline{\mathbf{n}}{\underline{\mathbf{n}}}^t-\epsilon d{\mathbf{n}}^{\times }-1,\varphi =\pi .$$

(49)

More generally, when extended to the dual setting, the composition law (2) yields

$$\langle {\underline{\mathbf{c}}}_2,{\underline{\mathbf{c}}}_1\rangle =\left(1+\epsilon {\lambda }_{\circ }\right)\langle {\mathbf{c}}_2,{\mathbf{c}}_1\rangle +\epsilon \lambda$$

(50)

with the notation

$${\displaystyle {\lambda }_{\circ }({\underline{\mathbf{c}}}_1,{\underline{\mathbf{c}}}_2)=\frac{{\mathbf{c}}_2·{\mathbf{d}}_1+{\mathbf{d}}_2·{\mathbf{c}}_1}{1-{\mathbf{c}}_2·{\mathbf{c}}_1},\lambda ({\underline{\mathbf{c}}}_1,{\underline{\mathbf{c}}}_2)=\frac{{\mathbf{d}}_2+{\mathbf{d}}_1+{\mathbf{c}}_2\times {\mathbf{d}}_1+{\mathbf{d}}_2\times {\mathbf{c}}_1}{1-{\mathbf{c}}_2·{\mathbf{c}}_1}}$$

and clearly, whenever ${\mathbf{c}}_2·{\mathbf{c}}_1=1$ one ends up with a transformation of the type (49).

If the condition ${\underline{r}}_{21}={\underline{g}}_{21}$ holds, we may use an alternative expression to (5) (see [5])

$${\underline{\tau }}_1={\displaystyle \frac{{\underline{\mathbf{a}}}_1\times {\underline{\mathbf{a}}}_2·\underline{\mathbf{c}}}{{\underline{g}}_{12}{\underline{\rho }}_1-{\underline{\rho }}_2}},{\underline{\tau }}_2={\displaystyle \frac{{\underline{\mathbf{a}}}_1\times {\underline{\mathbf{a}}}_2·\underline{\mathbf{c}}}{{\underline{g}}_{12}{\underline{\rho }}_2-{\underline{\rho }}_1}}$$

(51)

where ${\underline{\rho }}_i={\underline{\mathbf{a}}}_i·\underline{\mathbf{c}}$, to derive the solution more easily. These are further simplified in the case of orthogonal screw axes ${\underline{g}}_{12}=0$, while $Re\left({\underline{\rho }}_{1,2}\right)=0$ obviously leads to infinite elements, i.e., transformations of the type (49). Also, in the case ${\omega }_2\ne 0$ one has non-trivial closed rotation sequences (see [4]) which may be extended to closed rigid motion sequences with

$${\underline{\tau }}_1=\underline{\omega }{\underline{\Omega }}_{23}^{-1},{\underline{\tau }}_2=\underline{\omega }{\underline{\Omega }}_{31}^{-1},{\underline{\tau }}_3=\underline{\omega }{\underline{\Omega }}_{12}^{-1}$$

(52)

where ${\underline{\Omega }}_{ij}$ are the co-factors in the Gram determinant $\underline{\omega }={\underline{\omega }}_2$, e.g., ${\underline{\Omega }}_{23}={\underline{g}}_{12}{\underline{g}}_{31}-{\underline{g}}_{23}$ etc. This setting may be further generalized by allowing non-vanishing dual counterpart in the resulting transformation, i.e., factorizing a pure translation into three rigid motions with non-trivial rotational component. Neither of these is possible, however, in the degenerate case $\underline{\omega }=0$, although formally we have a solution for ${\underline{\mathbf{a}}}_3=\pm {\underline{\mathbf{a}}}_1$ with ${\underline{\tau }}_3=\mp {\underline{\tau }}_1$ and ${\underline{\tau }}_2=0$.

The gimbal lock singularity also persists in the dual setting, given by the condition

$${\underline{\mathbf{a}}}_3=\pm \underline{\mathcal{R}}{\underline{\mathbf{a}}}_1$$

(53)

in which case the first expression in (51) yields the dual scalar parameter for ${\underline{\varphi }}_1\pm {\underline{\varphi }}_3$, respectively. Let us consider a couple of examples before discussing the case of four factors.

Let us begin we a pair of screw axes given by ${\underline{\mathbf{a}}}_1={(0,1,-\epsilon )}^t$${\underline{\mathbf{a}}}_2={(1,0,\epsilon )}^t$, and a screw motion introduced with its dual Rodrigues’ vector representing a rotation1 by an angle $\varphi ={\displaystyle \frac{\pi }{3}}$ about a screw axis given by its direction vector $Re\left(\underline{\mathbf{c}}\right)$ and moment $Du\left(\underline{\mathbf{c}}\right)$. We calculate the matrix entries and find out that the condition ${\underline{r}}_{21}={\underline{g}}_{21}$ is satisfied, so a decomposition into a pair of screw motions exists with scalar parameters ${\tau }_{1,2}=1\pm \epsilon$ easily obtained using formula (51). Hence, one has $\underline{\mathbf{c}}=\langle {\underline{\mathbf{c}}}_2,{\underline{\mathbf{c}}}_1\rangle$ with ${\underline{\mathbf{c}}}_1={(0,1+\epsilon ,-\epsilon )}^t$ and ${\underline{\mathbf{c}}}_2={(1-\epsilon ,0,\epsilon )}^t$ as can easily be verified via (2), or in matrix terms, using (1) we can write

$$\mathcal{R}\left(\underline{\mathbf{c}}\right)=\left(\begin{array}{ccc}\hfill -2\epsilon & \hfill 0& \hfill 1\\ \hfill 1& \hfill 2\epsilon & \hfill 2\epsilon \\ \hfill -2\epsilon & \hfill 1& \hfill 0\end{array}\right)=\left(\begin{array}{ccc}\hfill 1& \hfill -\epsilon & \hfill \epsilon \\ \hfill \epsilon & \hfill \epsilon & \hfill -1\\ \hfill \epsilon & \hfill 1& \hfill \epsilon \end{array}\right)\left(\begin{array}{ccc}\hfill -\epsilon & \hfill \epsilon & \hfill 1\\ \hfill -\epsilon & \hfill 1& \hfill -\epsilon \\ \hfill -1& \hfill -\epsilon & \hfill -\epsilon \end{array}\right)=\mathcal{R}\left({\underline{\mathbf{c}}}_2\right)\mathcal{R}\left({\underline{\mathbf{c}}}_1\right).$$

Geometrically, ${\underline{\mathbf{c}}}_1$ and ${\underline{\mathbf{c}}}_2$ are interpreted as screw motions with quarter-turn rotational components and displacements $d_{1,2}=\pm 1$, as can be seen from (46). Next, we shall consider a decomposition into three factors $\mathcal{R}\left(\underline{\mathbf{c}}\right)=\mathcal{R}\left({\underline{\mathbf{c}}}_3\right)\mathcal{R}\left({\underline{\mathbf{c}}}_2\right)\mathcal{R}\left({\underline{\mathbf{c}}}_1\right)$ where ${\underline{\mathbf{a}}}_3={\underline{\mathbf{a}}}_1$ and $\mathcal{R}$ is rigid transformation of the type (49) with a screw line $\underline{\mathbf{n}}={\displaystyle \frac{1}{\sqrt{2}}}{(1+\epsilon ,1-\epsilon ,0)}^t$. Here the Davenport condition ${\underline{g}}_{12}={\underline{g}}_{13}=0$ is satisfied with $\underline{\Delta }=1$ and since ${\underline{\omega }}_i\equiv 0$ the solutions (3) are symmetric: ${\underline{\tau }}_1^{\pm }={\underline{\tau }}_3^{\pm }=\pm 1$ and ${\tau }_2^{\pm }=\pm (1+2\epsilon )$, which yields the matrix factorization

$$\left(\begin{array}{ccc}\hfill 2\epsilon & \hfill 1& \hfill 0\\ \hfill 1& \hfill -2\epsilon & \hfill 0\\ \hfill 0& \hfill 0& \hfill -1\end{array}\right)=\left(\begin{array}{ccc}\hfill 0& \hfill \epsilon & \hfill 1\\ \hfill -\epsilon & \hfill 1& \hfill -\epsilon \\ \hfill -1& \hfill -\epsilon & \hfill 0\end{array}\right)\left(\begin{array}{ccc}\hfill 1& \hfill -\epsilon & \hfill \epsilon \\ \hfill \epsilon & \hfill -2\epsilon & \hfill -1\\ \hfill \epsilon & \hfill 1& \hfill -2\epsilon \end{array}\right)\left(\begin{array}{ccc}\hfill 0& \hfill \epsilon & \hfill 1\\ \hfill -\epsilon & \hfill 1& \hfill -\epsilon \\ \hfill -1& \hfill -\epsilon & \hfill 0\end{array}\right)$$

along with its transpose (note that the matrix on the left is symmetric since $d=0$). Geometrically, the outer factor represents a quarter-turn about a screw axis given by its Plücker coordinates $Re\left({\underline{\mathbf{a}}}_1\right)$ and $Du\left({\underline{\mathbf{a}}}_1\right)$, while the one in the middle has a non-trivial displacement $d_2^{\pm }=\pm 2$ beside the quarter-turn rotation ${\varphi }_2^{\pm }=\pm {\displaystyle \frac{\pi }{2}}$, and its moment vector is $\mathbf{m}={(0,0,1)}^t$.

The transfer principle certainly allows us to do much more with dual matrices and quaternions. For interesting examples and possible applications we refer the reader to [9,10,11,25], and focus on generalizing what we referred to as the ’shift parameter’ construction.

4.2. Optimization of Screw Motions

We begin with the observation that formula (4) which provides the necessary and sufficient condition for factorization in $\mathrm{SO}\left(3\right)$ remains valid also in the case of generic rigid (screw) motions, since only the real part of the determinant actually remains under the square root as formula (44) shows. Hence, all inequalities concerning the relative angles, such as (11) or (28), also apply in this broader setting. It may seem counter-intuitive since pure translations cannot always be factorized unless the axes form a system of maximal rank, but that case goes beyond the validity of the transfer principle we rely on here, and yet, it is still possible in some cases, accounting for mutually cancelling rotational components in the decomposition. This can be considered as a variation of a closed rotation sequence solution which exists in two cases: non-coplanar axes, where formula (52) applies, and ${\underline{\mathbf{a}}}_3=\pm {\underline{\mathbf{a}}}_1$ with direct cancellation ${\tau }_3=\mp {\tau }_1$, ${\tau }_2=0$, leaving room for a residual slide, e.g.,

$$\langle {\mathbf{e}}_1+\epsilon {\mathbf{e}}_2,-{\mathbf{e}}_1+\epsilon {\mathbf{e}}_2\rangle =\epsilon \left({\mathbf{e}}_2+{\mathbf{e}}_3\right).$$

Having pointed that out, let us return to the regular case of screw motions with non-vanishing rotational component $\varphi \ne 0$. All equalities derived in the previous section generalize for the corresponding dual entities ${\underline{g}}_{ij}$, ${\underline{r}}_{ij}$, etc., while the inequalities apply for their real counterparts. Consider for example the decomposition (33) for a screw motion determined by its dual Rodrigues’ vector $\underline{\mathbf{c}}={(1,1+\epsilon ,1-\epsilon )}^t$ which is obviously a rotation by an angle $\varphi ={\displaystyle \frac{\pi }{3}}$ about a screw axis with Plücker coordinates $Re\left(\underline{\mathbf{c}}\right)$ and $Du\left(\underline{\mathbf{c}}\right)$. Let ${\tau }_i\left(\alpha \right)={\tau }_i^{+}$, so that ${\tau }_1\left(0\right)=1+\epsilon$, ${\tau }_2\left(0\right)={\displaystyle \frac{\epsilon }{2}}$ and ${\tau }_3\left(0\right)=1$, i.e., the decomposition involves a quarter-turn, a translation and a screw motion with ${\sum }_{i=1}^3\left|{\varphi }_i\right|=\pi$, ${\sum }_{i=1}^3\left|d_i\right|={\displaystyle \frac{3}{2}}$, namely

$$\left(\begin{array}{ccc}\hfill 0& \hfill \epsilon & \hfill 1\\ \hfill 1& \hfill \epsilon & \hfill 0\\ \hfill -\epsilon & \hfill 1& \hfill -\epsilon \end{array}\right)=\left(\begin{array}{ccc}\hfill 0& \hfill -1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\end{array}\right)\left(\begin{array}{ccc}\hfill 1& \hfill 0& \hfill \epsilon \\ \hfill 0& \hfill 1& \hfill 0\\ \hfill -\epsilon & \hfill 0& \hfill 1\end{array}\right)\left(\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill -\epsilon & \hfill -1\\ \hfill 0& \hfill 1& \hfill -\epsilon \end{array}\right)·$$

It is interesting to observe that in this case ${\tau }_1=1+\epsilon$ does not depend ot $\alpha$, while ${\tau }_2={\displaystyle \frac{\epsilon }{2cos\alpha }}$ and ${\tau }_3=1-\epsilon tan\alpha$, so the optimal factorization (w.r.t. both degrees of freedom) is at $\alpha =0$. Curiously enough, if we use the basis $\left\{{\underline{\mathbf{a}}}_i\right\}$ from our previous example, there is a stable factorization into a pair of screw motions $\underline{\mathbf{c}}=\langle (\epsilon -1){\underline{\mathbf{a}}}_3,(1+\epsilon ){\underline{\mathbf{a}}}_2\rangle$ utterly unaffected by $\alpha$.

Our next example illustrates how both the total rotation and total displacement in a screw decomposition can be affected independently by the additional parameter $\alpha$, which is kept real (angular) for the sake of representability. For our compound rigid motion we choose $\underline{\mathbf{c}}={\displaystyle \frac{1}{5}}{(3+4\epsilon ,4-3\epsilon ,\epsilon )}^t$ that is again a quarter-turn about a screw axis with direction vector $Re\left(\underline{\mathbf{c}}\right)$ and moment $Du\left(\underline{\mathbf{c}}\right)$. The decomposition takes place with respect to an orthonormal basis formed by ${\underline{\mathbf{a}}}_1$, ${\underline{\mathbf{a}}}_2$ and ${\underline{\mathbf{a}}}_3={\underline{\mathbf{a}}}_1\times {\underline{\mathbf{a}}}_2$ used in the previous example, so the dual extension of formula (33) applies here, and it is straightforward to assess the cost functions for the rotation and displacement components in the factorization separately (Figure 3). We should take into account that the minimum with respect to these two counterparts does not typically occur for the same value of $\alpha$, as can be seen in that example, so one has to prioritize according to the actual energy cost associated with these types of motion. This can also be done with the use ot statistical weights, adjusted for the situation.

Another observation worth mentioning is that while the intervals for $\alpha$ are preserved in the dual extensions, the solutions may have singularities in the displacement degree of freedom within these intervals. This should not be surprising, since the angular solutions also allow singularities: both due to the gimbal lock issue and in the form of infinite scalar parameters, corresponding to half-turns. Moreover, in the above example we impose an artificial restriction on $\alpha$ demanding that it is real, which breaks the rules of the game in a way. Typically, in the numerical experiments the overall displacement ’blows up’ when $\alpha$ is far from zero, but like with the angular sum, its minimum does not have to be at the origin.

4.3. A Brief Note on the Order of Finite Generation

Based on (10), (19) and (29) one may infer that the vertices of any (possibly degenerate) spherical triangle on ${\mathbb{S}}^2$ with perimeter at least $\pi$ can be used as axes in a decomposition into four factors for an arbitrary element of $\mathrm{SO}\left(3\right)$. In [24] this result is generalized as follows: if $\overline{\gamma }$ is the mean spherical distance of a broken closed geodesic on ${\mathbb{S}}^2$ connecting the vertices $\left\{{\mathbf{a}}_i\right\}$, then we have for the minimal number of rotations with invariant axes $\left\{{\mathbf{a}}_i\right\}$ two separate estimates given below (for an arbitrary axis and an arbitrary angle, respectively)

$$N_{\overline{\gamma }}\left(\varphi \right)\le 1+⌈{\displaystyle \frac{\left|\varphi \right|}{\overline{\gamma }}}⌉,N_{\overline{\gamma }}\left(\beta \right)\le 1+⌈{\displaystyle \frac{2\beta }{\overline{\gamma }}}⌉$$

(54)

where $\varphi$ is the angle of the compound rotation, $\beta$ is the spherical distance of its invariant axis to the closest vertex, and $\lceil ·\rceil$ the ’ceiling’ function (rounding from above). In particular, if the ${\mathbf{a}}_i$’s constitute a regular spherical polygon or a ’bouncing’ orbit between two axes (${\gamma }_{ij}=\overline{\gamma }$) and we want to cover the whole group (thus setting $\left|\varphi \right|=\pi$) the above reduces to an old result due to Lowenthal [23] on the order of finite generation of the group in $\mathrm{SO}\left(3\right)$.

Once more, since the discriminant condition (4) involves only the real (rotational) part of the screw motion, its direct consequences (10), (19) and (29) should also refer to the angular components. And again, one needs to be careful with the limitations of the transfer principle, discussed above, as to avoid falling into the common trap of wishful thinking.

5. Discussion

This paper proposes an closed for analytical solution to the optimization problem of rigid motion in attitude control of robots, drones and spacecraft demanding spatial rotations and translations, based on screw interactions modeled with projective dual quaternions (dual Rodrigues’ vectors). The goal is to minimize the energy output, measured by the total amount of motion in a sequence of four factors so that we have only one free parameter. The the Euclidean and the angular displacement come with weight coefficients in the cost function, since the difference in energy output in the two types of motion depends on the specific construction of the mechanism and environment in which it operates. The same refers also to rotations about different gimbals, which may be associated with different moments of inertia or external sources of asymmetry, resulting again in the introduction of weights in the cost function.

Note also that we do not consider here other optimization criteria, such as finding the smoothest trajectory interpolating a rigid motion for example, although this problem is also quite naturally approached using dual quaternions. Instead, we focus only on energy efficiency with just one additional parameter. This approach does not guarantee to provide the global minimum in all possible settings, but it is purely analytic, quite easy to implement, and still rather efficient: our calculations show it may reduce the cost of a maneuver by $30-40\%$ without the need of invoking complex numerical or machine learning algorithms used in other optimization methods (see [14,25,26]). Moreover, our solutions work with a relaxed version of the classical Davenport factorization condition for orthogonal axes, thus allowing for a wider variety of gimbal settings, e.g. tetrahedral or asymmetric, that might be useful in engineering applications.

Since our construction naturally extends to the spin covering group $\mathrm{SU}\left(2\right)$, it may be applied to qubits in quantum computing [16], computer vision and virtual reality [6]. A similar approach can be used for the proper Lorentz group in the plane ${\mathrm{SO}}^{+}(2,1)$ and its spin covering $\mathrm{SL}(2,\mathbb{R})\cong \mathrm{SU}(1,1)$ using the decomposition results obtained in [4]. In that setting, however, one needs to take into account, apart from the gimbal lock, also the isotropic singularity associated with the light cone in ${\mathbb{R}}^{2,1}$ which (unlike the former) is not affected by $\alpha$, being an invariant subspace independent on the the ${\tau }_i$’s. The ’dualization’ procedure works similarly in the hyperbolic setting, extending pseudo-rotations to general Poincaré motions in Minkowski $2+1$ space, but even within the homogeneous component, one has different types of transformations: pure rotations, Lorentz boosts, and isotropic displacements, which makes it harder to construct a meaningful cost function only from first principles, as factorizations in general may involve all three, e.g., the famous Iwasawa decomposition. It is an interesting direction in which this study can be generalized, although the $2+1$ Lorentz group has fewer applications compared to $\mathrm{SO}\left(3\right)$ or ${\mathrm{SO}}^{+}(3,1)$: apart from being used as a toy model for relativity, it plays a role in quantum mechanical scattering, the description of graphenes and certainly classical hyperbolic geometry. The entire Lorentz group, however, is a bit too wild to tame this way: although the factorizations may be described similarly using complex projective biquaternions, dualization does not have the same meaning and besides, the Plücker embedding is not trivial, which means we no longer deal with just plane transformations. However, interpreting the generalized ’axes’ $\left\{{\mathbf{a}}_k\right\}$ as generators of Abelian subgroups, the factorization problem still makes perfect sense in ${\mathrm{SO}}^{+}(3,1)\cong \mathrm{SO}(3,\mathbb{C})$ as well, and as long as one declares a meaningful cost optimization goal, a similar approach can also be applied to the (homogeneous) relativistic setting—from purely algebraic perspective, it is not so different to the one considered here (see [27,28]).