Are quaternions good or evil?

W.S. Harwin

race video

14/4/2021

$$\def\RR{{\bf R}} \def\bold#1{{\bf my #1}} \renewcommand{\vec}[1]{\underline{{#1}}}$$

# Good or evil?

## The story of Broom bridge, Dublin

Here as he walked by on the 16th of October 1843 Sir William Rowan Hamilton in a flash of genius discovered the fundamental formula for quaternion multiplication $i^2 = j^2 = k^2 = ijk = -1$ and cut it on a stone of this bridge.

• Quaternions were used by James Clerk Maxwell to describe EM theory; in 20 equations.
• Oliver Heaviside is responsible for the 4 equations known today as 'Maxwell's equations'
•  I came later to see that, as far as the vector analysis I required was concerned, the quaternion was not only not required, but was a positive evil of no inconsiderable magnitude; and that by its avoidance the establishment of vector analysis was made quite simple and its working also simplified, and that it could be conveniently harmonised with ordinary Cartesian work. ' Oliver Heaviside (1893)
•  Quaternions came from Hamilton after his really good work had been done; and, though beautifully ingenious, have been an unmixed evil to those who have touched them in any way, including Clerk Maxwell. ' W. Thompson, Lord Kelvin (1892)

# Coordinate transforms

• Needed for robotics, virtual reality, computer graphics etc.
• Each link can have a local coordinate frame
• Vector quantities (forces, positions, velocities and accelerations) can all be considered in a local coordinate frame
• Unit quaternions are ideal for considering vector rotations.
• Coordinate transforms can then separate into translation (vector addition) and rotation (quaternion multiplication)
• Homogeneous transforms
• Can use a single function (matrix multiply) for both translations and rotations
• Possible also with dual quaternions but there is no obvious benefit

# Vector rotations

If we have a vector $v$ we can introduce an operator $L_{\theta,\vec{u}}(v)$ rotates $v$ around axis defined by a unit vector $\vec{u}$ by an angle $\theta$.

This operator can be interpreted in two ways.

• Rotation of a vector by an angle $\theta$
• The same vector measured with respect to a second coordinate frame that has been rotated by $-\theta$

A rotational transform can be considered as

${}^b\vec{v}=L_{\theta,\vec{u}}({}^r\vec{v})$

a vector rotation of $\theta$ about $\vec{u}$, where ${}^b\vec{v}$ is the blue vector and ${}^r\vec{v}$ is the red vector (left).

Or it can be considered as

${}^n\vec{v}={}^n_sL_{-\theta,\vec{u}}({}^s\vec{v})$

that provides the coordinates of the vector in the new frame of reference rotated by $-\theta$ about the $z$ axis (right).

If $L_{-\theta,\vec{u}}$ has an inverse $L^{-1}_{-\theta,\vec{u}}=L_{\theta,\vec{u}}$ then

${}^s\vec{v}={}^n_sL^{-1}_{-\theta,\vec{u}}({}^n\vec{v})={}^s_nL_{\theta,\vec{u}}({}^n\vec{v})$

# Skew matrix arithmetic

If we have a vector

$\vec{u} = \begin{bmatrix} \omega _{x} & \omega _{y} & \omega _{z} \end{bmatrix}^T$

We can use a skew operator to form a skew matrix of the form

$S(\vec{u})=\tilde{S}=S_u= \begin{bmatrix} 0 & -\omega _{z} & \omega _{y}\\ \omega _{z} & 0 & -\omega _{x}\\ -\omega _{y} & \omega _{x} & 0 \end{bmatrix}$

$\tilde{S}$ or $S_u$ has the property that implements the cross product, i.e. $\vec{u}\times\vec{v}=\tilde{S}_u\vec{v}$ and $\vec{v}\times\vec{u}=-\tilde{S}_u\vec{v}=\tilde{S}_u^T\vec{v}$

It also has the following properties

\begin{align} S^T=-S\\ S^2-\vec{u}\vec{u}^T+I\vec{u}^T\vec{u}=0\\ S^3+S=(\vec{u}^T\vec{u}-1)S^T=(1-\vec{u}^T\vec{u})S \end{align}

If $\vec{u}$ is a unit vector then these properties simplify to

\begin{align} S^2-\vec{u}\vec{u}^T+I=0\label{eq:S4}\\ S^3=-S \label{eq:S5} \end{align}

Another property if $\vec{u}$ is a unit vector (from equation 4 and 5) of a unit vector skew matrix is that $S\vec{u}\vec{u}^T=0$

Further properties are

• S is singular if it has an odd dimension (Which is sad since it would make it easier to deal with exponential forms)
• (S+I) is invertable
• odd dimension skew matrices have only one real Eigenvalue which is 0.

# Rotation matrices and the exponential form

A property of $S_u$ is that the exponential will generate an orthogonal rotational matrix, usually expressed as

$$R=e^{S(\vec{u})\theta}=e^{S_{\vec{u}}\theta} \label{eq:expform}$$

where $\vec{u}$ is a unit vector and $\theta$ represents some rotation about $\vec{u}$.

This can be demonstrated by expanding the exponential as

$R=I + S\theta + S^2\frac{\theta^2}{2!} + S^3\frac{\theta^3}{3!} + S^4\frac{\theta^4}{4!} +...$$R=I + S\theta + S^2\frac{\theta^2}{2!} - S\frac{\theta^3}{3!} - S^2\frac{\theta^4}{4!} +...$$R=I + \left(S\theta - S\frac{\theta^3}{3!} +...\right) + \left(S^2\frac{\theta^2}{2!} - S^2\frac{\theta^4}{4!} +...\right)$$R=I + \sin\theta S + (1-\cos\theta)S^2$

This is sometimes called Euler's formula, possibly because it is an extension of $e^{i\theta}=\cos\theta+i\sin\theta$ but it is more properly called Rodrigues rotation formula.

Rodrigues formula was originally expressed in a vector notation, but using skew matrices is a more elegant approach.

# Some variants of Rodrigues rotation formula

Using the properties of skew matrices above the following variants of Rodrigues rotation formula are possible.

\begin{align} R = e^{S_u\theta} &= I + \sin\theta S_u + (1-\cos\theta)S_u^2\\ &= I + \sin\theta S_u + S_u^2 -\cos\theta S_u^2\\ &= I + \sin\theta S_u+(1-\cos\theta)(\vec{u}\,\vec{u}^T - I)\\ &= I\cos\theta + \sin\theta S_u + (1-\cos\theta)\vec{u}\,\vec{u}^T\\ &= \sin\theta S_u + \vec{u}\,\vec{u}^T + (I-\vec{u}\,\vec{u}^T)\cos\theta \end{align}

and evidently

$R^T= e^{-S_u\theta}= e^{S^T_u\theta} =I-\sin\theta S_u + (1-\cos\theta) S_u^2$

# Differentiation of a rotation matrix

if $S_u$ is based on a unit vector and constant and $R=e^{S_u\theta}$ (equation 6) then

$\frac{dR}{d\theta}=S_u e^{S_u\theta} = S_u R$

and

$\dot{R}=\frac{dR}{dt}=\frac{dR}{d\theta}\frac{d\theta}{dt}=S_u \left(R\frac{d\theta}{dt}\right)=\omega S_u R$

where $\omega=\frac{d\theta}{dt}$

Often the Skew matrix is considered as representing an angular velocity with components in a particular coordinate frame, that is $\Omega=S_u \omega$ represents rotation about a unit vector $\vec{u}$ at an angular speed of $\omega$.

# Quaternions vs rotation matrices

Quaternions (when used to calculate rotations) are a 4-tuple, i.e. a set of 4 Real numbers.

• Efficient to calculate
• Good for applying splines and Bezier curves (only 4 parameters to spline)
• Fewer terms to consider when doing calculations
• Easy to interchange with matrices (e.g. via Rodrigues rotation formula)
• Would manage the Dirac belt trick/plate trick ( SU(2) (e.g. quaternions) double-covers SO(3) (the rotations group))
• Extends to affine transforms with dual quaternions

• Efficient notation ${}^av= {}^a_bR\,\,{}^bv$ against ${}^av= {}^a_bq\,\, {}^bv\,\, {}^a_bq^{-1}$
• More widely used in publications, easier to understand
• Possible to interchange with angle axis and quaternions. (axis of the rotation comes via the Eigen vectors, the angle of rotation comes via the matrix trace)
• extends to affine transforms with homogeneous transforms.

• Easy to conceptualise
• Easy to interchange with matrices and quaternions

# Quaternions

## Notation abuse

• Quaternion multiplication is non commutative and is sometimes explicitly represented by $l\otimes q$ and sometimes not.
• Binary operations (multiplication and addition) will be assumed if one of the variables is a quaternion (e.g. a vector pretending to be a quaternion $(0,\vec{v})$, or a real number $(4,0,0,0)$ )
• We will ignore the $ijk$ multiplications and use the scalar/vector representation
• Thus for the tuple $r=(r_0,r_1,r_2,r_3)$ we will represent as $r=(r_0,\vec{r}$) with $\vec{r}=[r_1\ r_2\ r_3]^T$ (note this may annoy some mathematicians)
• Occasionally we may also abuse the $+$ and $-$ signs and write quaternions as $r=r_0+\vec{r}$ this is because quaternions can be written as $r=r_0+r_1i+r_2j+r_3k$.

# Quaternion basics

## Quaternion conjugate

$q^* = q_0 -\vec{q_v}= (q_0, -\vec{q_v})$

## Quaternion norm

$\|q\| = q_0^2+ \vec{q_v}\cdot \vec{q_v}$

## Quaternion multiplication

Given

$r = (r_0,\vec{r_v})$$q = (q_0, \vec{q_v})$

we get

\begin{align*} q \otimes r = & (q_0 + \vec{q}_v)( r_0 + \vec{r}_v) = (q_0,\vec{q}_v)\otimes ( r_0, \vec{r}_v) \\ = &(q_0r_0 - \vec{q}_v{\bf.}\vec{r}_v), (q_0\vec{r}_v + r_0\vec{q}_v + \vec{q}_v\times\vec{r}_v) \end{align*}

where $\times$ is the vector cross product.

## Quaternion inverse

$q^{-1} = \frac{q_0 -\vec{q}_v}{q_0^2+ \vec{q}_v\cdot\vec{q}_v}$

which can be expressed as the conjugate over the squared norm.

## unit quaternions

Any quaternion with a norm of 1. A useful unit quaternion is

$q=\left(\cos\left(\frac\theta{2}\right),\vec{u}\sin\left(\frac\theta{2}\right)\right)$

where $u$ is a unit vector.

# Quaternion rotation

Rotation of a quaternion representing a vector $v$ about a unit vector $\vec{u}$ through an angle $(\theta)$ is given by $$q' = rqr^{-1}$$ where

$r=\cos(\theta/2) + \sin(\theta/2) \vec{u}$

and the quaternion $q=(0,\vec{v})$ is a vector subset of the quaternion group.

Expanding the generalised rotation in vector notation (and using the conjugate rather than the inverse to saved work)

$rqr^* = q_0(r_0^2+\vec{r_v}{\bf .}\vec{r_v}) \,,\, (r_0^2 - \vec{r_v}{\bf .}\vec{r_v})\vec{q_u}+2r_0(\vec{r_v}\times\vec{q_u}) +2(\vec{q_u}{\bf .}\vec{r_v})\vec{r_v}$

if $q_0=0$

$rqr^* = (r_0^2 - \vec{r_v}{\bf .}\vec{r_v})\vec{q_u}+2r_0(\vec{r_v}\times\vec{q_u}) +2(\vec{q_u}{\bf .}\vec{r_v})\vec{r_v}$

which is a vector

This can be rewritten (expanding the half angles to full angles) as

$$rqr^* = \cos(\theta)\vec{q_u} + \sin(\theta)(\vec{r_v}\times\vec{q_u})+(1-\cos(\theta))(\vec{q_u}{\bf .}\vec{r_v})\vec{r_v} \label{eq:rod99}$$

rewriting the quaternion elements

$\displaystyle L_{\theta,\vec{k}}(\vec{v}) = \cos(\theta)\vec{v} + \sin(\theta)(\vec{k}\times\vec{v})+(1-\cos(\theta))(\vec{v}{\bf .}\vec{k})\vec{k}$

The equation 12 is the vector form of Rodrigues rotation formula

source

# Quaternions as a matrix

$q\otimes r$ can be done as matrix product $Q_L r$ if $Q_L$ is a 'left product' matrix and $r$ is considered as a 4 element vector.

The conjugate is simply the matrix transpose so $q^*\,r$ can be done as $Q_L^T r$

$Lmatrix= Q_L= \begin{bmatrix} q_{0} & -q_{1} & -q_{2} & -q_{3}\\ q_{1} & q_{0} & -q_{3} & q_{2}\\ q_{2} & q_{3} & q_{0} & -q_{1}\\ q_{3} & -q_{2} & q_{1} & q_{0} \end{bmatrix} = \begin{bmatrix} q_0 & -\vec{q}_v^T\\ \vec{q}_v & S_{\vec{q}}+I q_0 \end{bmatrix}$

A `right product' matrix is needed and $r\otimes q$ is then done as the calculation $r^T Q_R$ where

$Rmatrix= Q_R= \begin{bmatrix} q_0 & \vec{q}_v^T\\ -\vec{q}_v & S_{\vec{q}}+I q_0 \end{bmatrix}$

Note that $Q^T_L\ne Q_R$

It is then possible to calculate

$r'=qrq*=Q_L(r^T Q_R^T)^T$

The resulting 4x4 matrix $Q_L Q_R$ embeds the 3x3 rotation matrix in the lower right corner.

# Quaternion differentiation

Consider a quaternion changing over an incremental period of time from $q(t)$ to $q(t+\Delta t)$. This incremental change is given by

$q(t+\Delta t) =q(\Delta t)\,\, q(t)$

If $q(t+\Delta t)$ and $q(t)$ are unit quaternions then $q(\Delta t)$ must be a unit quaternion. Assume that in the incremental time $\Delta t$ the angle of $q(\Delta t)$ changes by some $\Delta\theta$

so define the quaternion differential to be

$\frac{d q}{d t}=\dot{q}=\lim_{\Delta t\to 0}\frac{q(t+\Delta t)-q(t)}{\Delta t}\\ =\lim_{\Delta t \to 0}\frac{q(\Delta t)\,\,q(t)-q(t)}{\Delta t}\\ =\lim_{\Delta t \to 0}\frac{(q(\Delta t)\,-\,1)q(t)}{\Delta t}$

If we assume $q(\Delta t)$ is of the form

$q(\Delta t)=\left(\cos\left(\Delta\theta/2\right)+\vec{u}\sin\left(\Delta\theta/2\right)\right)$

Then as $\Delta t$ approaches 0, $q(\Delta\theta)$ approaches $(1+\vec{u}\,\Delta\theta/2)$ and hence

$\frac{d q}{d t} =\lim_{\Delta \to 0}\frac{(1+\vec{u}\Delta\theta/2\,-\,1)q(t)}{\Delta t} =\lim_{\Delta \to 0}\frac{\vec{u}\Delta\theta/2}{\Delta t}q(t) =\frac12\vec{u}\frac{d\theta}{d t}q(t)$

If $q$ defines a vector or coordinate frame then the differential can be considered as an angular velocity around unit $\vec{u}$ at a speed of $\frac{d\theta}{d t}$ so

$\frac{d q}{d t} =\frac12 \vec{\omega} q(t)$

and $\vec\omega$ is the angular velocity vector.

A more complete explanation given in [Jia2020]

# So good or evil?

Vote now.

## Good

• Quaternions use 4 parameters to specify a rotation, vs 9 for rotation matrices, or 4 when expressed in exponential form.
• Easily described as a (half) angle and a rotation axis
• Numerically easy to normalise when compared with rotation matrices, or exponential form.
• important to keep integration errors under control
• Easy to change between quaternions, rotation matrices, and angle axis.
• Fewer calculations to compute a rotation
• perhaps not important given overheads of operating systems and speeds of GPUs and CPUs.

## Evil

• Longer to write out, rotations need 3 symbols vs 2.
• Are not suitable for every problem.
• e.g. recursive dynamics well established using rotation matrices.
• Don't have a cool exponential form

## Conclusion.

• Probably neither.
• But not a universal panacea.
• In effect, a good compliment to rotation matrices and angle-axis representations.
• Don't overlook Rodrigues forms

# Refs

Shuster, M. 1993, "A Survey of Attitude Representations", Journal of the Astronautical Sciences, 41(4):349-517 See equations and discussion in the paper above, p463-464. (http://www.ladispe.polito.it/corsi/Meccatronica/02JHCOR/2011-12/Slides/Shuster_Pub_1993h_J_Repsurv_scan.pdf)

Yan-Bin Jia, Problem Solving Techniques for Applied Computer Science, Quaternions Lecture course Sep 3, 2020 (https://web.cs.iastate.edu/~cs577/) (https://web.cs.iastate.edu/~cs577/handouts/quaternion.pdf)

Quaternions, Interpolation and Animation Erik B. Dam Martin Koch Martin Lillholm (https://web.mit.edu/2.998/www/QuaternionReport1.pdf)

Lee, Byung-Uk (1991). Unit Quaternion Representation of Rotation - Appendix A, Differentiation with Quaternions - Appendix B (PDF) (Ph. D. Thesis). Stanford University (http://home.ewha.ac.kr/~bulee/quaternion.pdf)

## Bibs

Plenty of details at Wolfram Mathworld, Wikipedia, Numberphile (although surprisingly quiet about differentiation)

Possibly a good explanation! (https://fgiesen.wordpress.com/2012/08/24/quaternion-differentiation/) google:Differentiation with respect to the rotation quaternion

# Appendix

## Differential of Rodrigues form

$dR=(\cos\theta S_u +\sin\theta S_u^2)d\theta$