# Why is SLERP not the same as this method for averaging quaternions?

I understand that the weighted average of at least 2 quaternions is as described here in equation 13. To summarise, the steps are:

1. 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".
2. 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]

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 + np.sin((1 - weight) * theta) * qs) / 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)):.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")


Essentially, because equally spaced points on a chord of a circle don't subtend equal arcs on the circle. Here's an exaggerated version to illustrate the effect: The points along the chord at $$y=1/2$$ are equally spaced, but you can see that the angles they make aren't equal. Slerping between two quaternions corresponds to an angular interpolation, where the slerp parameter is the portion of the angle to the interpolated point; linear interpolation (and in this case the weighted sum) corresponds to the chord here, of exactly equal distances between the points, but not equal angles.

• I feel like this is an explanation for why I can't just straight up interpolate the quaternion coefficients. But I'm not doing that. In fact, the outputs of the method I described are unit quaternions. Am I missing your point though? Nov 11, 2022 at 19:23
• @AlexanderSoare I believe what you're doing effectively amounts to normalizing the interpolated results, or something very similar to that; the core point is that taking weighted averages is similar to taking linear sections along the chord between the two quaternion points. Nov 11, 2022 at 19:29
• @AlexanderSoare I believe in this case you're actually taking what amounts to a linear interpolation between matrix elements of SO(3), but the core point is the same; SO(3) isn't a 'flat' space, so linear interpolations between the Euclidean representation of its elements don't correspond to linear interpolation of 'distance' (angle) in the actual SO(3) metric space. Nov 11, 2022 at 19:30
• Hm it's not clear to me that computing the weighted sum of outer products of quaternions and then finding the eigenvector that corresponds to the largest eigenvalue of that, is the same is doing linear interpolation then normalizing the result. But then again, I'm not so sure how to interpret what it is that I'm doing. All I know is that it's the most cited method I've found for averaging quaternions. (this is in response to your second last message). Nov 11, 2022 at 19:31
• You may also be interested in mathweb.ucsd.edu/~sbuss/ResearchWeb/spheremean which gives a weighted average on the sphere (which includes the quaternion sphere, in this case) that they claim does a better job of respecting spherical distance. Nov 11, 2022 at 20:20

Consider the axis-angle representation of quaternions $$\mathbf{q}_1$$ and $$\mathbf{q}_2$$, \begin{aligned} \mathbf{q}_1 &= \cos\left(\frac{\theta_1}{2}\right) + \hat{n}_1 \sin\left(\frac{\theta_1}{2}\right) \\ \mathbf{q}_2 &= \cos\left(\frac{\theta_2}{2}\right) + \hat{n}_2 \sin\left(\frac{\theta_2}{2}\right) \\ \end{aligned} representing rotation by angle $$\theta_k$$ around unit axis vector $$\hat{n}_k$$.

Direct interpolation by $$\lambda$$ followed by normalization yields \begin{aligned} \mathbf{q}_u &= (1 - \lambda) \mathbf{q}_1 + \lambda \mathbf{q}_2 \\ \mathbf{q} &= \frac{\mathbf{q}_u}{\left\lVert\mathbf{q}_u\right\rVert} \\ \end{aligned}

Consider the case where the rotation axes are perpendicular, $$\hat{n}_1 \perp \hat{n}_2$$. The angle $$\varphi$$ that the rotation axis makes with $$\hat{n}_1$$ is $$\varphi = \arctan\left(\frac{\lambda}{1 - \lambda}\right)$$ so that at $$\lambda = 0$$, $$\varphi = 0°$$, and at $$\lambda = 1$$, $$\varphi = 90°$$ (since the two rotation axes are perpendicular to each other). This is not linear, but it is close: For the rotation part, the rotation angle $$\theta$$ is $$\theta = 2 \arccos\biggl( (1 - \lambda) \cos\left(\frac{\theta_1}{2}\right) + \lambda \cos\left(\frac{\theta_2}{2}\right) \biggr)$$ If you compare this to $$\theta = (1 - \lambda) \theta_1 + \lambda \theta_2$$, you'll find that this is not linear either, but is rather close. For $$t_1 = 0°$$, $$t_2 = 180°$$, which shows the maximum possible effect. The smaller $$\lvert \theta_2 - \theta_1 \rvert$$, the closer to linear you get.

This means that normalized weighted sum of two quaternions is a valid interpolation between the two quaternions, but it is only an approximation of the interpolation that minimizes the curve length that any point on the unit sphere will form, when interpolating from one quaternion to the other. It is surprisingly close, though. (It could be... I just don't think it is, exactly.)

In fact, I am not exactly convinced that SLERP provides that minimum either. What it does provide, is an axis that tracks the great circle that intersects with both axes, with the angular change linear in $$\lambda$$, plus a linear change in rotation angle $$\theta$$ with respect to $$\lambda$$. (In other words, the lines in black in the two above diagrams.)

You can argue which one is more correct, but you really need to define "correct" first. The direct interpolation produces results acceptable to human perception or space games, but for a spacecraft, you use a different approach, one depending on exactly what you want to minimize by your approach.

• Thanks for this. But I do mention specifically in my question that I am not asking about directly averaging quaternion coefficients. I am referring to a commonly used method for "averaging" quaternions and I outline this in my question. Nov 21, 2022 at 8:40