Many times it is necessary to calculate the weighted average of 2 rotors (I am restricting this discussion to 3D Euclidean rotations). That’s easy in Geometric Algebra, given two rotors and we can solve for knowing that for the equation holds (identity at zero) and for we have (the final rotor) so:
thus we have:
Which satisfy our requirements and is smooth. Finally we get the expression: for getting rotors from to .
Notice you don’t need to know the angle as in the SLERP formula, however it can be shown that is the same as the SLERP interpolation (you get the same result). So we can actually interpolate rotors using SLERP if convenient.
But what about interpolating 3, 4 or rotors? it turns out the interpolation of more than 2 rotors is not that easy, we can extend our previous approach, but that’s only one approach out many and there are trade-offs we need to be aware of. So let’s begin with the extension of the previous approach.
The Rotor Manifold
3D rotors live in a manifold i.e., a curved space, which is basically a four dimensional sphere called . Also the 3D rotors form a group called SO(3) of orthogonal isometries (i.e., transformations that preserve length and angles on the sphere and also preserve orientation). The SO(3) group is also a Lie group i.e., the rotor exponential is a map from tangent space to the manifold i.e., exponential map can project a tangent vector into the manifold, or in other words it can project an euclidean bivector to a rotor (so euclidean bivectors are tangent vectors and 3D rotors are points of the manifold). However the rotor exponential only works in such way at the north pole of i.e., at the identity. That means if we have a rotor which moves from the identity rotor to some “point” rotor in , the exponential works just fine. But if you have a rotor which moves from an arbitrary point rotor to other arbitrary point rotor then the exponential map cannot be used. Fortunately the Lie group SO(3) is “translation invariant” which means that we can translate an arbitrary point rotor to the north pole of , use the exponential map there and then translate back to point rotor . And the beauty of Lie groups is that all those translations can be done using rotors.
If we have rotors we can attempt to average the rotors in the “tangent space” at the north pole of .
But what does it mean? computation wise that amounts to do: where . We can see this is basically the average of bivectors (tangent vectors) in the tangent space at the north pole of . Then we are using the exponential map to project the averaged result to the manifold so we get an averaged rotor. Notice that we are assuming that all rotors have their tangent space at the north pole of i.e., that they are all “connected to the identity” by a geodesic curve (a curve departing from the identity rotor to the given rotor).
That’s a reasonable assumption for many applications. For 2 rotors it will give an interpolation close to SLERP although without constant velocity. It is also cheap to compute (closed form) and one can do high order rotor interpolation using polynomial basis functions, B-Spline or any other.
The geodesic distance covered by the geodesic curve in can be measured as which is basically . Note that I am using the logarithm defined at the north pole of i.e., at the identity. The quality of an average can be measured in terms of minimizing the distance i.e., the rotor average is the one at the “center” or at “equal distance” to all input rotors, so a sensible measure can be where $R$ is the average rotor and is our measure of the quality (the less is the better is the quality of the average).
We can improve the Log-Euclidean average by doing a gradient descent on . Fortunately there is a simple way to do that. The “Weiszfeld Algorithm” is a known optimization algorithm for finding the geometric mean on a Riemannian manifold. Note that here is a scalarfield. The gradient vector field of is , which is a great result (we derive it in next section). The gradient points in the negative direction of the geodesic path minimizing . The goal is to find a mean rotor such that . That minimization makes sense since the gradient belongs to the tangent space at the identity where all combinations are linear i.e., the of each point is translated to the identity i.e., . The gradient descent iteration is as follows:
- Set and repeat steps 2-4 until
As you can see in step 3, we project the gradient to the manifold and translate the point from the identity (north pole) to where it belongs moving along the geodesic curve. The rotor is the average rotor minimizing the geodesic distance. Instead of a gradient descent one can do a Newton type optimization, minimizing the squared distance . That improved version is the “Karcher Algorithm“.
Derivation of the L1 Gradient for Weiszfeld Algorithm
Having geodesic curve from to parameterized by the scalar given by and a function measuring the geodesic distance from to , where and the covariant derivative is given by:
by definition of
Replacing in the original equation:
Since we can identify any tangent vector along the geodesic with then the following inner product represent the differential :
By identification with we have
So the gradient of is and we can apply that to by linearity to get the gradient of it.
Time Derivative of
by definition of
By commutativity of exponentials of the same bivector.
After evaluating at
By linearization of
Chord Metric Minimizer Average
Other way to define the average is in terms of the Euclidean distance i.e., the average that minimize the so called “chordal distance” or distance of rotors embedded in space (4-dimensional vector space), defined as:
Unfortunately for rotors the distance is not well defined since a rotor and represent the same rotation and it doesn’t matter which one to choose, so the distance should be zero but it isn’t. So the mean rotor minimizing chordal distance on average is meaningful only if you choose all the to have the same “orientation”.
On the other hand, the chord metric for matrices under Frobenius norm is well defined and it can be shown, using Rodrigues formula, that chordal length is related to the sine of the angle as
For rotors you can also minimize the angle using the following distance function:
where the norm is the Euclidean norm in the embedding space of 3D bivectors . So this distance is measurimg the same quantity than the chord distance of matrices under Frobenius norm (up to a scalar). We will attempt to minimize the squared distance since we can a closed form expression. Basically we have to solve the following optimization problem:
This functional is minimizing the squared angle
Terminology: lets denote the rotor as the sum of the bivector and the scalar so . Similarly the rotor .
Where denote the commutator product of bivectors. This can be trivially expressed to matrix form by embedding bivectors in and changing commutator product by cross product of vectors.
Where is the cross product matrix and is the identity matrix.
Let’s denote the matrix made exclusively with terms of rotor , then we can express the quadratic form as:
The matrix has the following form:
The whole optimization problem can be expressed as:
The quadratic form is minimized when is equal to the eigenvector of associated with the min eigenvalue.
This problem we solved above can be restated as
which is maximizing . This dual problem convenient since this the quadratic form associated is much simpler than the primal problem:
Recall that and this is equivalent to:
Notice the negative sign of is absorbed by the Euclidean inner product so there is no mistake in the sign. So the quadratic form is:
where and so the quadratic form is:
The quadratic form is maximized when is equal to the eigenvector of associated with the max eigenvalue.
The above quadratic form was derived by Markley et al http://www.acsu.buffalo.edu/~johnc/ave_quat07.pdf in 2007 for finding quaternion’s average and it is quite popular in the graphics industry since it is simple to implement and fast to compute i.e., closed form. However Markley et al paper lacks clarity, they used a much more complex mathematical approach to derive the quadratic form, based on matrices related to the Wahba’s Problem, probably because the first author was familiar to those methods. Anyway, they also gave an obscure interpretation to the solution. They assert that this solution is a “maximum likelihood estimate” of an attitude matrix which is totally unrelated to the problem i.e., the Fisher information matrix of matrix errors in SO(3). I found that interpretation totally clueless and irrelevant. The real interpretation is that it minimizes the squared angle of rotors (or quaternions) embedded in the vector space, since it is equivalent to primal problem. You can also say it minimizes the chordal distance of the equivalent matrices under Frobenius norm.