Rotating a single vector using a quaternion

By February 9, 2019

This is a rehash of something I wrote in a forum post something like 10 years ago. It turns out that forum prohibited archiving in its robots.txt so it’s not on the Internet Archive. The original operator of said forum hosted a copy (with broken formatting) of the original posts for a while, but at a different URL, breaking all links. And now that copy’s gone too, again with archiving disabled apparently. Sigh.

I got asked about this yesterday; I don’t have a copy of my original derivation anymore either, but here’s my best attempt at reconstructing what I probably wrote back then, and hopefully it won’t get lost again this time.

Suppose you’re given a unit quaternion q and a vector v. Quaternions are a common rotation representation in several fields (including computer graphics and numerical rigid-body dynamics) for reasons beyond the scope of this post. To apply a rotation to a vector, one computes the quaternion product q v q^*, where v is implicitly identified with the quaternion with real (scalar) part 0 and v as its imaginary part, and q^* denotes the conjugate of q. Such quaternions with a real part of 0 are also referred to as “pure imaginary” quaternions. Anyway, the result of the above product is another pure imaginary quaternion, corresponding to the rotated vector.

This is all explained and motivated elsewhere, and I won’t bother doing so here. Now generally, you often want to apply the same transformation to many vectors. In that case, you’re better off turning the quaternion into the equivalent 3×3 rotation matrix first. You can look up the formula elsewhere or work it out from the expression above with some algebra. That’s not the topic of this post either.

But sometimes you really only want to transform a single vector with a quaternion, or have other reasons for not wanting to expand to an explicit 3×3 (or larger) matrix first, like for example trying to minimize live register count in a shader program. So let’s look at ways of doing that.

The direct way

First, I’m going to split quaternions in their real (scalar) and imaginary (vector) parts, writing q=(q_r, q_{xyz}) where q_r is the real part and q_{xyz} imaginary. The conjugate of q is given by q^*=(q_r, -q_{xyz}).

The product of two quaternions a and b is given by

ab = (a_r b_r - a_{xyz} \cdot b_{xyz}, a_r b_{xyz} + b_r a_{xyz} + a_{xyz} \times b_{xyz})

where \cdot denotes the usual dot product and \times the cross product. If you’re not used to seeing this in vector notation, I recommend using this (or something more abstract like geometric algebra) for derivations; writing everything in terms of scalars and the i, j, k basis elements makes you miss the forest for the trees.

Anyway, armed with this formula, we can compute our final product without too much trouble. Let’s start with the qv part:

qv = (-q_{xyz} \cdot v, q_r v + q_{xyz} \times v)

Not so bad. Now we have to multiply it from the right by q^*, which I’ll do in multiple steps. First, let’s take care of the real part, by multiplying our just-computed values for qv with q^* using the general multiplication formula above:

(qvq^*)_r = -q_r (q_{xyz} \cdot v) - ((q_r v + q_{xyz} \times v) \cdot (-q_{xyz}))
= -q_r (q_{xyz} \cdot v) + q_r (v \cdot q_{xyz}) + ((q_{xyz} \times v) \cdot q_{xyz})
= q_{xyz} \cdot (q_{xyz} \times v) = 0

because the first two dot products are identical and the cross product of q_{xyz} and v is perpendicular to q_{xyz}. This proves that qvq^* is indeed a pure imaginary quaternion again, just like the v we started out with.

Nice to know, but of course we’re actually here for the vector part:

(qvq^*)_{xyz} = (-q_{xyz} \cdot v) (-q_{xyz}) + q_r (q_r v + q_{xyz} \times v)
+ (q_r v + q_{xyz} \times v) \times (-q_{xyz})
= (q_{xyz} \cdot v) q_{xyz} + q_r^2 v + q_r (q_{xyz} \times v) + q_{xyz} \times (q_r v + q_{xyz} \times v)
= (q_{xyz} \cdot v) q_{xyz} + q_r^2 v + 2 q_r (q_{xyz} \times v) + q_{xyz} \times (q_{xyz} \times v)

which so far has used nothing fancier than the antisymmetry of the cross product a \times b = -b \times a.

If we pull out and name the shared cross product, we get

u = q_{xyz} \times v
(qvq^*)_{xyz} = (q_{xyz} \cdot v) q_{xyz} + q_r^2 v + 2 q_r u + q_{xyz} \times u

which is the direct expression for the transformed vector from the formula. This is what you get if you just plug everything into the formulas and apply basic algebraic simplifications until you run out of obvious things to try (which is exactly what we did).

In terms of scalar operation count, this boils down to two cross products at 6 multiplies and 3 additions (well, subtractions) each; one 3D dot product at 3 multiplies and 2 additions; 3 vector-by-scalar multiplications at 3 multiplies each; two scalar multiplies (to form q_r^2 and 2q_r, although the latter can also be computed via addition if you prefer); and finally 3 vector additions at 3 adds each. The total operation count is thus 26 scalar multiplies and 17 additions, unless I miscounted. For GPUs, a multiply-add “a*b+c” generally counts as a single operation, and in that case the scalar operation count is 9 scalar multiplies and 17 scalar multiply-adds.

Stand back, I’m going to try algebra!

We can do better than that. So far, we haven’t used that q is a unit quaternion, meaning that q_r^2 + q_{xyz} \cdot q_{xyz} = 1. We can thus plug in (1 - q_{xyz} \cdot q_{xyz}) for q_r^2, which yields:

(qvq^*)_{xyz} = (q_{xyz} \cdot v) q_{xyz} + (1 - q_{xyz} \cdot q_{xyz}) v + 2 q_r u + q_{xyz} \times u
= v + (q_{xyz} \cdot v) q_{xyz} - (q_{xyz} \cdot q_{xyz}) v + 2 q_r u + q_{xyz} \times u

This might look worse, but it’s progress, because we can now apply the vector triple product identity

a \times (b \times c) = (a \cdot c)b - (a \cdot b)c

with a=q_{xyz}, b=q_{xyz}, c=v, leaving us with:

(qvq^*)_{xyz} = v + q_{xyz} \times (q_{xyz} \times v) + 2 q_r u + q_{xyz} \times u
= v + q_{xyz} \times u + 2 q_r u + q_{xyz} \times u
= v + 2 q_r u + 2 q_{xyz} \times u

The two remaining terms involving u both multiply by two, so we use a slightly different shared temporary for the final version of our formula:

t = 2 q_{xyz} \times v
(qvq^*)_{xyz} = v + q_r t + q_{xyz} \times t

Final scalar operation count: without multiply-adds, the two cross products sum to 12 multiplies and 6 adds, scaling t by two takes 3 multiplies (or adds, your choice), and the final vector-by-scalar multiply and summation take 3 multiplies and 6 adds, respectively. That’s 18 multiplies and 12 adds total, or 15 multiplies and 15 adds if you did the doubling using addition.

Either way, that’s a significant reduction on both counts.

With multiply-adds, I count a total of 3 multiplies and 9 multiply-adds for the cross products (the first cross product has nothing to sum to, but the second does); either 3 multiplies or 3 adds (your choice again) for the doubling; and another 3 multiply-adds for the v + q_r t portion. That’s either 6 multiplies and 12 multiply-adds or 3 multiplies, 3 adds and 12 multiply-adds. Furthermore some GPUs can fold a doubling straight into the operand without counting as another operation at all; in that case we get 3 multiplies and 12 multiply-adds.

For comparison, multiplying a vector by a 3×3 rotation matrix takes 9 multiplies and 6 additions (or 3 multiplies plus 6 multiply-adds). So even though this is much better than we started out with, it’s generally still worthwhile to form that rotation matrix explicitly if you plan on transforming lots of vectors by the same quaternion, and aren’t worried about GPU vector register counts or similar.