Using biquaternions for Lorentz transformations

What’s the problem?

A long time ago, while debugging my game, I found a curious issue when using single-precision floating point numbers (floats) rather than double-precision floating point numbers (doubles). The game crashed on various assertions (checks I have in the game to make sure it is doing the “right thing”.)

These crashes were due to numerical errors accumulating step after step during the rigid body calculations. This issue is similar to the one arising when dividing 1 by 3. The result is 0.3333333… However, you cannot store infinites 3-s inside a computer. You have to truncate the number at some point. The truncation introduces small errors that can accumulate and grow.

Well, matrices (I use them to represent Lorentz transformations in the game) have even a worse problem: the set of all possible matrices is typically too big for what they are needed for. For example, rotations can be represented using a particular  kind of 3 by 3 matrices. That’s pretty cool as composing two rotations ends up corresponding to multiplying the two respective rotation matrices. Unfortunately, this is true only in the ideal world of maths, where multiplications are not affected by rounding errors. In the real world, however, the product of two rotation matrices is not guaranteed to be a rotation matrix. This can become a problem if many matrices are multiplied to obtain a final one.

Unfortunately, this is what I was doing: I was using 4 by 4 matrices, one for each rigid body, to represent Lorentz transformations. At each step I was multiplying the  “current” matrix by another Lorentz matrix. The numerical error was thus growing with time.

Why not solve the issue by using doubles?

You may argue that the solution to this problem is doing the calculations with double precision. Well, this is how I dealt with the issue until now. But this is not a real fix. It’s a matter of time until the error grows and becomes big enough to cause troubles. It also feels uncomfortable having numerical errors that accumulate and have no recipe to “cure” them.

Why quaternions and biquaternions?

Here is where quaternions can help. While a rotation matrix uses 3x3 = 9 numbers to represent a rotation, a quaternion only uses 4. And, actually, 4 is still too many. A rotation only needs 3 numbers to be uniquely determined. However, a generic quaternion is easier to correct back to a “rotation quaternion”: dividing all its components by its norm is all that it takes. This “correction operation” can then be applied periodically, to ensure the quaternion doesn’t “deteriorate” with time.

What about Lorentz transformations? Well, it turns out that Lorentz transformations are kind-of a generalization of rotations. It isn’t surprising, then, that an appropriate generalization of quaternions can be used to do something similar. Biquaternions are indeed this generalization. They are like quaternions, except that they are made of four complex numbers rather than four real numbers.

How did it go?

I started the work by refactoring the C++ object that I used to store 4x4 matrices, to make it easier to reimplement it using biquaternions. (e.g. I removed all explicit accesses to the matrix entries) I then wrote a new implementation of it using biquaternions. I used this paper as a reference.

Finally, I played the old trick of temporarily changing the code so that it would use both implementations and cross check them at each step. When doable, this gives you pretty good confidence that the new implementation is correct.

That’s it. I can now switch between the two implementations at compile time.

The positives…

The switch to biquaternions allows me to keep numerical errors under control. I can now use single-precision floating-point numbers for rigid body calculations.

Also, Lorentz transformations require less memory (8 floats, rather than 16 doubles - a 4x size reduction.)

Composition between Lorentz transformations is also considerably more efficient, as it only requires one biquaternion-multiplication.

The negatives…

However, applying a Lorentz transformation to a four-vector is faster using matrices. The biquaternion representation does indeed require two biquaternion multiplications, which is considerably slower than doing a simple matrix-vector multiplication. This is unfortunate as these operations are more frequently carried out than compositions.

A solution to this problem would be to use both representations, i.e. convert a biquaternion to a matrix whenever there is a need to Lorentz-transform multiple four-vectors. I’d like to implement this change, but… not now. The next step for me is to focus on the demo release. Further work on the biquaternion support can wait, unless I find it can't :)