Invert 4x4 matrix - Numerical most stable solution

Question

I want to invert a 4x4 matrix. My numbers are stored in fixed-point format (1.15.16 to be exact).

With floating-point arithmetic I usually just build the adjoint matrix and divide by the determinant (e.g. brute force the solution). That worked for me so far, but when dealing with fixed point numbers I get an unacceptable precision loss due to all of the multiplications used.

Note: In fixed point arithmetic I always throw away some of the least significant bits of immediate results.

So - What's the most numerical stable way to invert a matrix? I don't mind much about the performance, but simply going to floating-point would be to slow on my target architecture.

Answer 1

You might consider doubling to 1.31 before doing your normal algorithm. It'll double the number of multiplications, but you're doing a matrix invert and anything you do is going to be pretty tied to the multiplier in your processor.

For anyone interested in finding the equations for a 4x4 invert, you can use a symbolic math package to resolve them for you. The TI-89 will do it even, although it'll take several minutes.

If you give us an idea of what the matrix invert does for you, and how it fits in with the rest of your processing we might be able to suggest alternatives.

-Adam

Answer 2

Meta-answer: Is it really a general 4x4 matrix? If your matrix has a special form, then there are direct formulas for inverting that would be fast and keep your operation count down.

For example, if it's a standard homogenous coordinate transform from graphics, like:

[ux vx wx tx]
[uy vy wy ty]
[uz vz wz tz]
[ 0  0  0  1]

(assuming a composition of rotation, scale, translation matrices)

then there's an easily-derivable direct formula, which is

[ux uy uz -dot(u,t)]
[vx vy vz -dot(v,t)]
[wx wy wz -dot(w,t)]
[ 0  0  0     1    ]

(ASCII matrices stolen from the linked page.)

You probably can't beat that for loss of precision in fixed point.

If your matrix comes from some domain where you know it has more structure, then there's likely to be an easy answer.

Answer 3

To minimize truncation errors and other badness, use "pivoting" - see the chapter on inverting matrices in Numerical Recipes. They have the best explanation i've found so far.

Answer 4

If the matrix represents an affine transformation (many times this is the case with 4x4 matrices so long as you don't introduce a scaling component) the inverse is simply the transpose of the upper 3x3 rotation part with the last column negated. Obviously if you require a generalized solution then looking into Gaussian elimination is probably the easiest.

Answer 5

I think the answer to this depends on the exact form of the matrix. A standard decomposition method (LU, QR, Cholesky etc.) with pivoting (an essential) is fairly good on fixed point, especially for a small 4x4 matrix. See the book 'Numerical Recipes' by Press et al. for a description of these methods.

This paper gives some useful algorithms, but is behind a paywall unfortunately. They recommend a (pivoted) Cholesky decomposition with some additional features too complicated to list here.

Answer 6

Plain old Gaussian elimination would work well.

It depends on what libraries/classes/structures you're using. You could take a look at the GSL.