I understand that the weighted average of at least 2 quaternions is as described here in equation 13. To summarise, the steps are:
- Take the weighted sum $\mathbf{M} = \sum_i w_i \mathbf{q}_i \mathbf{q}_i^\intercal$, where $\mathbf{q}$ is the vector of quaternion coefficients, $i$ indexes the set of quaternions and their corresponding weights $w_i$. In words: "take the weighted sum of outer products of quaternions with themselves".
- Find the eigenvector corresponding to the greatest eigenvalue of $\mathbf{M}$. That is the result.
To compare that against SLERP we restrict ourselves to only two quaternions and have weights $\{\alpha, 1 - \alpha \}$.
Being a programmer it was faster for me to write code to compare that with SLERP (see appendix for Python code), and I found that the two methods give the same answer if $\alpha \in \{0, 0.5, 1.0\}$ but otherwise not. In fact, SLERP gives the "right" answer intuitively, whereas the weighted average method does not. Would someone please point me in the direction of an explanation as to why the weighted average does not seem to be doing what one would naively think it should? ie the same as SLERP.
Edit: To be clear, when I say "averaging quaternions", I am never referring to straight linear averaging, and it's clear to me why the latter wouldn't work. I am only referring to the method described above and referenced from the link.
Edit 2: Probably the best clue I have is equation 17 in the link. It shows that the algorithm for averaging quaternions is equivalent to minimizing the weighted sum of the squares of the sines of the half-error-angles between the inputs and the dummy quaternion.
On the other hand, I suppose SLERP effectively minimizes the weighted sum of error-angles squared (I could be wrong, please correct me if so). Or written with the same notation as the paper we'd be computing:
$$ \text{argmin}_{\mathbf{q}} w_i \delta \phi^2 $$
I suppose then that at this point, my question distills down to: why would anyone call equation 13 of the linked paper a "weighted average" when that's not what it seems to be doing intuitively? And why is it okay to use this method to do things like state estimation by taking weighted averages of measurements?
Appendix: Python code comparing weighted average of quaternions and SLERP
from typing import Tuple, Sequence
import random
import numpy as np
from scipy.spatial.transform import Rotation
from scipy.spatial.transform import Slerp
def qwxyz_to_single_cover(q_wxyz: np.ndarray) -> np.ndarray:
""" Maps quaternions to a subspace that represents a single-cover on SO(3).
Args:
q_wxyz: (4,) array representing a scalar-first quaternion.
Returns:
(4,) array for the remapped quaternion.
"""
for c in q_wxyz:
if c < 0:
q_wxyz = -q_wxyz
break
elif c > 0:
break
return q_wxyz
def quaternion_weighted_sum(Q: Sequence[np.ndarray], weights: Sequence[float]) -> np.ndarray:
""" Quaternion weighted mean as per http://www.acsu.buffalo.edu/~johnc/ave_quat07.pdf.
Args:
Q: a (n, 4) array of scalar-first quaternions
weights: a (n,) array of weights. They need not sum to 1.
Return:
A (4,) array for the resulting quaternion.
"""
A = sum(w * (np.outer(q, q)) for q, w in zip(Q, weights))
# Get the eigenvector corresponding to largest eigenvalue.
return np.linalg.eigh(A)[1][:, -1]
def slerp(qs: Tuple[np.ndarray, np.ndarray], weight: float):
""" SLERP
Args:
qs: A tuple of two (4,) arrays representing scalar-first quaternions
weights: A weight for the interpolation. The interpolation looks like
weight * q1 + (1 - weight) * q2.
Returns:
A (4,) array for the resulting quaternion.
"""
theta = np.arccos(np.dot(q1, q2))
return (np.sin(weight * theta) * qs[0] + np.sin((1 - weight) * theta) * qs[1]) / np.sin(theta)
def axis_angle_from_quaternion(q: np.ndarray) -> Tuple[np.ndarray, float]:
""" Given a quaternion get the axis of rotation and angle of rotation in radians.
"""
r = Rotation.from_quat(np.roll(q, -1)).as_rotvec()
norm = np.linalg.norm(r)
# Force angle to be in [-pi, pi)
if norm >= np.pi:
norm = 2 * np.pi - norm
return r / norm, norm
def print_q_info(q, postfix):
ax, ang = axis_angle_from_quaternion(q)
print(f"ax{postfix}: {ax}, ang{postfix}: {np.rad2deg(ang):.1f} deg, q{postfix}: {q}")
for i in range(3):
### Make two quaternions.
# q1 = qwxyz_to_single_cover(np.roll(Rotation.random().as_quat(), shift=1))
# First quaternion is trivial.
q1 = np.array([1, 0, 0, 0])
# Second quaternion is random.
q2 = qwxyz_to_single_cover(np.roll(Rotation.random().as_quat(), shift=1))
print_q_info(q1, 1)
print_q_info(q2, 2)
### Choose a weight for interpolation.
weight = 0.75 #random.random()
print(f"weight: {weight}")
# We expect the angle of rotation to be a linear interpolation.
print(f"expected ang: {(1 - weight) * np.rad2deg(axis_angle_from_quaternion(q2)[1]):.1f} deg")
print()
### Compute and print.
print("Quaternion weighted sum")
q_ws = quaternion_weighted_sum((q1, q2), (weight, 1 - weight))
print_q_info(q_ws, "_ws")
print()
print("Slerp")
q_sl = slerp((q1, q2), weight)
print_q_info(q_sl, "_sl")
# print("Scipy slerp")
# sl = Slerp([0, 1], Rotation.from_quat(np.roll(np.vstack([q1, q2]), shift=-1, axis=-1)))
# q_SL = np.roll(sl(1 - weight).as_quat(), shift=1)
# print_q_info(q_SL, "_SL")
print("\n====\n")