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

Question

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][:, -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")

Answer 1

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.

Answer 2

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.