Implementing proper normalization for unit dual quaternions is not straightforward. I got it wrong the first time and I got it wrong the second time I tried it, but this tought me more about dual numbers and dual quaternions, which I want to summarize in this article. The correct normalization procedure can be very simple in the end. It is a case in which several pages of calculations lead to a concise and simple solution that requires about 4 lines of code. Most of the things discussed here were extracted from this pull request.
Background¶
Let's first establish some fundamentals. A dual number $r + d \epsilon$ consists of a real part $r$ and a dual part $d$ with $\epsilon^2 = 0$ and $\epsilon \neq 0$. A quaternion $r + \boldsymbol{v}$ consists of a real part $r \in \mathbb{R}$ and a vector part $\boldsymbol{v} \in \mathbb{R}^3$. A dual quaternion is a combination of the two. A dual quaternion $\sigma = p + q \epsilon$ consists of a real quaternion $p = r_p + \boldsymbol{v}_p$ and a dual quaternion $q = r_q + \boldsymbol{v}_q$. Normalizing a dual quaternion means the dual quaternion should have unit norm, that is, $\sigma \cdot \sigma^* = 1$ (with the dual quaternion product $\cdot$). It is often assumed that it is enough to ensure that the real quaternion has unit norm by dividing the whole dual quaternion by the norm of the real quaternion. This works for most cases, but it is an approximation. We see this immediately, when we write the unit norm property correctly: $\sigma \cdot \sigma^* = 1 + 0 \epsilon$. Ensuring that the dual part is 0 is not trivial. There is a great solution that I found on stackoverflow here, but it lacks some details that I would like to write down here.
The following definitions and properties of quaternions, dual quaternions, and dual numbers are required to understand the solution.
Quaternion Product¶
$$q_1 \cdot q_2 = r_1 r_2 - \boldsymbol{v}_1^T \boldsymbol{v}_2 + r_1 \boldsymbol{v}_2 + r_2 \boldsymbol{v}_1 + \boldsymbol{v}_1 \times \boldsymbol{v}_2$$Quaternion Conjugate¶
$$q^* = (r + \boldsymbol{v})^* = r - \boldsymbol{v}$$Quaternion Product with Conjugate¶
From the last two definitions, we can see easily that
$$q \cdot q^* = r^2 + \boldsymbol{v}^T\boldsymbol{v} - r \boldsymbol{v} + r \boldsymbol{v} - \boldsymbol{v} \times \boldsymbol{v} = r^2 + \boldsymbol{v}^T \boldsymbol{v},$$which is just a real number without the vector part.
Dual Quaternion Product¶
$$\sigma_1 \cdot \sigma_2 = (p_1 + q_1 \epsilon) \cdot (p_2 + q_2 \epsilon) = p_1 \cdot p_2 + (p_1 \cdot q_2 + p_2 \cdot q_1) \epsilon$$Dual Quaternion Conjugate¶
There are three different dual quaternion conjugates. We will use the quaternion conjugate of the two individual quaternions, which is the inverse in the case of unit dual quaternions.
$$\sigma^* = (p + q \epsilon)^* = p^* + q^* \epsilon$$Norm of a Dual Quaternion¶
Now we can start computing the (squared) norm of a dual quaternion. We need the dual quaternion product and the dual quaternion conjugate:
$$\sigma \cdot \sigma^* = (p + q \epsilon) \cdot (p + q \epsilon)^* = p \cdot p^* + (p q^* + q p^*) \epsilon$$We can further break this down with $p = r_p + \boldsymbol{v}_p$ and $q = r_q + \boldsymbol{v}_q$ to
$$p \cdot p^* = r_p^2 + \boldsymbol{v}_p^T \boldsymbol{v}_p$$$$p q^* = r_p r_q + \boldsymbol{v}_p^T \boldsymbol{v}_q - r_p \boldsymbol{v}_q + r_q \boldsymbol{v}_p - \boldsymbol{v}_p \times \boldsymbol{v}_q$$$$q p^* = r_q r_p + \boldsymbol{v}_q^T \boldsymbol{v}_p - r_q \boldsymbol{v}_p + r_p \boldsymbol{v}_q - \boldsymbol{v}_q \times \boldsymbol{v}_p$$With $\boldsymbol{v}_p \times \boldsymbol{v}_q = - \boldsymbol{v}_q \times \boldsymbol{v}_p$, we can furthermore simplify (underlined terms cancel each other) $$p q^* + q p^* = (r_p r_q + \boldsymbol{v}_p^T \boldsymbol{v}_q \underline{- r_p \boldsymbol{v}_q + r_q \boldsymbol{v}_p} - \underline{\boldsymbol{v}_p \times \boldsymbol{v}_q}) + (r_q r_p + \boldsymbol{v}_q^T \boldsymbol{v}_p \underline{- r_q \boldsymbol{v}_p + r_p \boldsymbol{v}_q} - \underline{\boldsymbol{v}_q \times \boldsymbol{v}_p}) = 2 (r_p r_q + \boldsymbol{v}_p^T \boldsymbol{v}_q),$$
which does not contain a vector part anymore.
To summarize, the squared norm of a dual quaternion is a dual number of the form
$$\sigma \cdot \sigma^* = r_p^2 + \boldsymbol{v}_p^T \boldsymbol{v}_p + 2 (r_p r_q + \boldsymbol{v}_p^T \boldsymbol{v}_q) \epsilon,$$that is, the vector parts of the quaternions are $0$.
When we write $p = (r_p, v_{1p} v_{2p}, v_{3p})^T$ and $q = (r_q, v_{1q} v_{2q}, v_{3q})^T$ as 4d vectors, $r_p^2 + \boldsymbol{v}_p^T \boldsymbol{v}_p = p^T p = ||p||_2^2$ is the squared Euclidean norm of the real quaternion and $r_p r_q + \boldsymbol{v}_p^T \boldsymbol{v}_q = p^T q$ is the dot product of the real and dual quaternion, hence
$$\sigma \cdot \sigma^* = p^T p + 2 p^T q \epsilon.$$Example¶
The following function will compute the squared norm. It returns an array of shape (2,) with the real part in the first component and the dual part in the second component.
import pytransform3d.transformations as pt
def dual_quaternion_squared_norm(sigma):
return pt.concatenate_dual_quaternions(
sigma, pt.dq_q_conj(sigma, unit=False), unit=False
)[[0, 4]]
We can also use our simplified expression to implement the norm:
def dual_quaternion_squared_norm_simplified(sigma):
p = sigma[:4]
q = sigma[4:]
real = p @ p
dual = 2 * p @ q
return np.array([real, dual])
Let's try to generate an arbitrary 8d vector, interpret it as a dual quaternion and compute its norm.
import numpy as np
rng = np.random.default_rng(42)
sigma = rng.normal(size=8)
sigma
dual_quaternion_squared_norm(sigma)
dual_quaternion_squared_norm_simplified(sigma)
Normalization Procedure (Direct Solution)¶
To normalize a dual quaternion, we multiply it with the inverse square root of its squared norm.
Inverse Square Root of a Dual Number¶
Computing the inverse square root of a dual number is not trivial. One derivation is outlined here. However, I find the following more easy to understand.
We define the inverse square root of a dual number $s + t \epsilon$ as
$$\frac{1}{\sqrt{s + t \epsilon}} = u + v \epsilon$$and square the expression on both sides to
$$\frac{1}{s + t \epsilon} = u^2 + 2uv \epsilon.$$We multiply by $s + t \epsilon$ and get
$$1 + 0 \epsilon = (u^2 + 2uv \epsilon)(s + t \epsilon) = u^2s + (u^2t + 2uvs) \epsilon,$$which we can split into
$$1 = u^2s \Leftrightarrow u = \pm\frac{1}{\sqrt{s}}$$$$0 \epsilon = (u^2t + 2uvs) \epsilon \Leftrightarrow v \epsilon = - \frac{0.5 t}{(\pm \sqrt{s})^3} \epsilon$$Normalization¶
Hence, to normalize a dual quaternion, we multiply the dual quaternion with the dual number
$$\frac{1}{\sqrt{s + t \epsilon}} = \pm\left(\frac{1}{\sqrt{s}} - \frac{0.5 t}{(\sqrt{s})^3} \epsilon\right).$$We will only use the positive solution here, multiplying this normalization factor with the dual quaternion gives us the solution
$$(p + q \epsilon) \left(\frac{1}{\sqrt{s}} - \frac{0.5 t}{(\sqrt{s})^3} \epsilon\right) = \frac{p}{\sqrt{s}} + \left(\frac{q}{\sqrt{s}} - \frac{0.5 t}{(\sqrt{s})^3} p \right) \epsilon.$$Plugging in $s = p^T p$ and $t = 2 p^T q$ from the dual quaternion norm, we would get the normalization factor
$$\frac{1}{\sqrt{s + t \epsilon}} = \pm\left(\frac{1}{\sqrt{p^T p}} - \frac{p^Tq}{\left(\sqrt{p^Tp}\right)^3} \epsilon\right),$$hence,
$$(p + q \epsilon) \left(\frac{1}{\sqrt{s}} - \frac{0.5 t}{(\sqrt{s})^3} \epsilon\right) = \frac{p}{\sqrt{p^T p}} + \left(\frac{q}{\sqrt{p^T p}} - \frac{p^Tq}{\left(\sqrt{p^Tp}\right)^3} p \right) \epsilon.$$Side Note: Two Solutions¶
Why are there two solutions? As you can see, the difference is the sign of the normalization factor. Multiplying any unit dual quaternion by $-1$ gives us a unit dual quaternion that represents exactly the same transformation just like it is the case for unit quaternions and rotations.
Example¶
def norm_dual_quaternion(sigma):
sigma_prod = pt.concatenate_dual_quaternions(
sigma, pt.dq_q_conj(sigma, unit=False), unit=False
)
s = sigma_prod[0]
t = sigma_prod[4]
u = 1.0 / np.sqrt(s)
v = -0.5 * t / np.sqrt(s)**3
p = sigma[:4] * u
q = sigma[4:] * u + sigma[:4] * v
return np.hstack((p, q))
Again, we can write down a simplified version that we derived above.
def norm_dual_quaternion_simplified(sigma):
p = sigma[:4]
q = sigma[4:]
s = p @ p
u = 1.0 / np.sqrt(s)
v = -(p @ q) / np.sqrt(s)**3
p = sigma[:4] * u
q = sigma[4:] * u + sigma[:4] * v
return np.hstack((p, q))
sigma_unit = norm_dual_quaternion(sigma)
sigma_unit
dual_quaternion_squared_norm(sigma_unit).round(3)
sigma_unit = norm_dual_quaternion_simplified(sigma)
sigma_unit
dual_quaternion_squared_norm(sigma_unit).round(3)
We see that both functions convert an arbitrary 8d vector to a unit dual quaternion with squared norm $1 + 0 \epsilon$.
Edge Case¶
It could also happen that the real quaternion $p = 0$. In this case, $u$ is not defined. We will omit this case here, but the solution would be to set $p = 1$ and then proceed with the normalization procedure.
norm_dual_quaternion(np.zeros(8))
Normalization in Two Steps¶
To simplify the terms involved in the normalization procedure, we can split the normalization up into two steps:
- Multiply $\sigma$ with the normalization factor $\frac{1}{\sqrt{s}} = \frac{1}{\sqrt{p^T p}}$ for the real quaternion to obtain $\sigma' = \frac{\sigma}{\sqrt{s}}$, which ensures that the norm is $1 + t' \epsilon$.
- Multiply $\sigma' = p' + q' \epsilon$ with the normalization factor $(1 - \frac{t'}{2} \epsilon)$ for the dual quaternion (since s' = 1) to obtain $\sigma'' = p' + q' \epsilon - p'^T q' p' \epsilon$.
def norm_dual_quaternion_two_steps(sigma):
p = np.copy(sigma[:4])
q = np.copy(sigma[4:])
# 1. ensure unit norm of the real quaternion
real_norm = np.sqrt(p @ p)
p /= real_norm
q /= real_norm
# 2. ensure orthogonality
q -= p @ q * p
return np.hstack((p, q))
sigma_unit = norm_dual_quaternion_two_steps(sigma)
sigma_unit
dual_quaternion_squared_norm(sigma_unit).round(3)
The solution is simple. The way to get there was not so simple.
Geometrical Interpretation¶
Considering that unit dual quaternion $$\sigma = p + q \epsilon = p + 0.5 (0 + \boldsymbol{t}) p \epsilon,$$ represents a rigid transformation with a rotation quaternion $p$ and a translation vector $\boldsymbol{t}$, the procedure makes sense geometrically. We have to divide both quaternions by the norm of the real quaternion because it is part of both quaternions $p$ and $q$. In the second step, we remove the part of $q$ that is not orthogonal to $p$: the vector projection of $q$ on $p$, which is $p^T q p$ when $p^T p = 1$.
Mistakes¶
So what were the errors that I saw / made when I implemented the normalization procedure? Here is a list of the most important ones.
- Ignoring the orthogonality constraint. In the beginning I thought it would be enough to divide the dual quaternion by the norm of the real quaternion. This is even suggested in several tutorials, such as this. I guess it is good enough as an approximation for most use cases, but sometimes it is not.
- Dividing only the real quaternion by its norm and forgetting about the dual quaternion. Both quaternions should be devided by the norm of the real quaternion.
- Forgetting to take the square root of the real component $s$ of the squared dual quaternion norm. Since we never compute the actual norm $\sqrt{\sigma \sigma^*}$, but only the squared norm $\sigma \sigma^*$ of a dual quaternion $\sigma$ in the direct solution, this has to be kept in mind.
- Using the updated real quaternion to update the dual quaternion in the direct solution. This is necessary in the two-step solution, but it leads to confusing errors in the direct solution.