39

I was hoping someone can point out an efficient formula for 4x4 affine matrix transform. Currently my code uses cofactor expansion and it allocates a temporary array for each cofactor. It's easy to read, but it's slower than it should be.

Note, this isn't homework and I know how to work it out manually using 4x4 co-factor expansion, it's just a pain and not really an interesting problem for me. Also I've googled and came up with a few sites that give you the formula already (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm). However this one could probably be optimized further by pre-computing some of the products. I'm sure someone came up with the "best" formula for this at one point or another?

3
  • 1
    Why not use some existing libraries? Chances are those are already optimized.
    – kennytm
    Apr 12, 2010 at 18:37
  • 5
    That is true. Unfortunately that matrix code is in Java and then compiled by GWT. Most libraries simply won't work. Also it's a fairly narrow application. I'm just dealing with 4x4 matrices. I don't want to link a huge linear algebra library just to get inverse() and multiply() functionality.
    – Budric
    Apr 12, 2010 at 18:54
  • 1
    @kennytm what if we want to learn? Apr 9, 2022 at 17:32

5 Answers 5

54

You should be able to exploit the fact that the matrix is affine to speed things up over a full inverse. Namely, if your matrix looks like this

A = [ M   b  ]
    [ 0   1  ]

where A is 4x4, M is 3x3, b is 3x1, and the bottom row is (0,0,0,1), then

inv(A) = [ inv(M)   -inv(M) * b ]
         [   0            1     ]

Depending on your situation, it may be faster to compute the result of inv(A) * x instead of actually forming inv(A). In that case, things simplify to

inv(A) * [x] = [ inv(M) * (x - b) ]
         [1] = [        1         ] 

where x is a 3x1 vector (usually a 3D point).

Lastly, if M represents a rotation (i.e. its columns are orthonormal), then you can use the fact that inv(M) = transpose(M). Then computing the inverse of A is just a matter of subtracting the translation component, and multiplying by the transpose of the 3x3 part.

Note that whether or not the matrix is orthonormal is something that you should know from the analysis of the problem. Checking it during runtime would be fairly expensive; although you might want to do it in debug builds to check that your assumptions hold.

Hope all of that is clear...

5
  • So the first formula you got from "blockwise inversion" (en.wikipedia.org/wiki/Invertible_matrix)? Or is there another mental trick? I'm not quite sure about the inv(A) * x = inv(M) * (x - b). First, they're different sizes - do you remove a row from A on the left or add a row on the right? Second, I'm not sure how that equation comes about. Third, I'm not sure what you're solving for in that equation. Oliver keeps mentioning not computing inverses symbolically, but I don't know what that means - I need the inverse to do inverse transform. If you have time I'd like to know.
    – Budric
    Apr 13, 2010 at 15:32
  • 1
    I edited the inv(A) * x formula to make the dimensions clearer. The first formula was from en.wikipedia.org/wiki/Affine_transformation. Forgetting about affine transforms for a minute, in general, when you're solving Ax = b, you want the solution inv(A)*b. But often times you don't need to actual form inv(A), just compute what the product *would be. Back to affine transforms, in 3D applications, you might not actually need the inverse of the matrix, you just want the inverse transform acting on (multiplying) a vector. If that's the case, then using the formula might be faster.
    – celion
    Apr 13, 2010 at 18:22
  • 1
    Even if you do need to store the matrix inverse, you can use the fact that it's affine to reduce the work computing the inverse, since you only need to invert a 3x3 matrix instead of 4x4. And if you know that it's a rotation, computing the transpose is much faster than computing the inverse, and in this case, they're equivalent.
    – celion
    Apr 13, 2010 at 18:24
  • And here's a better explanation of what I meant by computing inv(A) * x: johndcook.com/blog/2010/01/19/dont-invert-that-matrix (the author is a frequent poster on SO).
    – celion
    Apr 13, 2010 at 18:38
  • Clear now. Very interesting.
    – Budric
    Apr 13, 2010 at 21:54
22

To follow-up on pkhaler's and Robin Hilliard's excellent responses above, here is Robin's ActionScript 3 code converted into a C# method. Hopefully this can save some typing for other C# developers, as well as C/C++ and Java developers in need of a 4x4 matrix inversion function:

public static double[,] GetInverse(double[,] a)
{
    var s0 = a[0, 0] * a[1, 1] - a[1, 0] * a[0, 1];
    var s1 = a[0, 0] * a[1, 2] - a[1, 0] * a[0, 2];
    var s2 = a[0, 0] * a[1, 3] - a[1, 0] * a[0, 3];
    var s3 = a[0, 1] * a[1, 2] - a[1, 1] * a[0, 2];
    var s4 = a[0, 1] * a[1, 3] - a[1, 1] * a[0, 3];
    var s5 = a[0, 2] * a[1, 3] - a[1, 2] * a[0, 3];

    var c5 = a[2, 2] * a[3, 3] - a[3, 2] * a[2, 3];
    var c4 = a[2, 1] * a[3, 3] - a[3, 1] * a[2, 3];
    var c3 = a[2, 1] * a[3, 2] - a[3, 1] * a[2, 2];
    var c2 = a[2, 0] * a[3, 3] - a[3, 0] * a[2, 3];
    var c1 = a[2, 0] * a[3, 2] - a[3, 0] * a[2, 2];
    var c0 = a[2, 0] * a[3, 1] - a[3, 0] * a[2, 1];

    // Should check for 0 determinant
    var invdet = 1.0 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    var b = new double[4, 4];

    b[0, 0] = ( a[1, 1] * c5 - a[1, 2] * c4 + a[1, 3] * c3) * invdet;
    b[0, 1] = (-a[0, 1] * c5 + a[0, 2] * c4 - a[0, 3] * c3) * invdet;
    b[0, 2] = ( a[3, 1] * s5 - a[3, 2] * s4 + a[3, 3] * s3) * invdet;
    b[0, 3] = (-a[2, 1] * s5 + a[2, 2] * s4 - a[2, 3] * s3) * invdet;

    b[1, 0] = (-a[1, 0] * c5 + a[1, 2] * c2 - a[1, 3] * c1) * invdet;
    b[1, 1] = ( a[0, 0] * c5 - a[0, 2] * c2 + a[0, 3] * c1) * invdet;
    b[1, 2] = (-a[3, 0] * s5 + a[3, 2] * s2 - a[3, 3] * s1) * invdet;
    b[1, 3] = ( a[2, 0] * s5 - a[2, 2] * s2 + a[2, 3] * s1) * invdet;

    b[2, 0] = ( a[1, 0] * c4 - a[1, 1] * c2 + a[1, 3] * c0) * invdet;
    b[2, 1] = (-a[0, 0] * c4 + a[0, 1] * c2 - a[0, 3] * c0) * invdet;
    b[2, 2] = ( a[3, 0] * s4 - a[3, 1] * s2 + a[3, 3] * s0) * invdet;
    b[2, 3] = (-a[2, 0] * s4 + a[2, 1] * s2 - a[2, 3] * s0) * invdet;

    b[3, 0] = (-a[1, 0] * c3 + a[1, 1] * c1 - a[1, 2] * c0) * invdet;
    b[3, 1] = ( a[0, 0] * c3 - a[0, 1] * c1 + a[0, 2] * c0) * invdet;
    b[3, 2] = (-a[3, 0] * s3 + a[3, 1] * s1 - a[3, 2] * s0) * invdet;
    b[3, 3] = ( a[2, 0] * s3 - a[2, 1] * s1 + a[2, 2] * s0) * invdet;

    return b;
}
0
22

Just in case someone would like to save some typing, here's an AS3 version I wrote based on page 9 (more efficient version of Laplace Expansion Theorem) of the link posted above by phkahler:

public function invert() : Matrix4 {
    var m : Matrix4 = new Matrix4();

    var s0 : Number = i00 * i11 - i10 * i01;
    var s1 : Number = i00 * i12 - i10 * i02;
    var s2 : Number = i00 * i13 - i10 * i03;
    var s3 : Number = i01 * i12 - i11 * i02;
    var s4 : Number = i01 * i13 - i11 * i03;
    var s5 : Number = i02 * i13 - i12 * i03;

    var c5 : Number = i22 * i33 - i32 * i23;
    var c4 : Number = i21 * i33 - i31 * i23;
    var c3 : Number = i21 * i32 - i31 * i22;
    var c2 : Number = i20 * i33 - i30 * i23;
    var c1 : Number = i20 * i32 - i30 * i22;
    var c0 : Number = i20 * i31 - i30 * i21;

    // Should check for 0 determinant

    var invdet : Number = 1 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    m.i00 = (i11 * c5 - i12 * c4 + i13 * c3) * invdet;
    m.i01 = (-i01 * c5 + i02 * c4 - i03 * c3) * invdet;
    m.i02 = (i31 * s5 - i32 * s4 + i33 * s3) * invdet;
    m.i03 = (-i21 * s5 + i22 * s4 - i23 * s3) * invdet;

    m.i10 = (-i10 * c5 + i12 * c2 - i13 * c1) * invdet;
    m.i11 = (i00 * c5 - i02 * c2 + i03 * c1) * invdet;
    m.i12 = (-i30 * s5 + i32 * s2 - i33 * s1) * invdet;
    m.i13 = (i20 * s5 - i22 * s2 + i23 * s1) * invdet;

    m.i20 = (i10 * c4 - i11 * c2 + i13 * c0) * invdet;
    m.i21 = (-i00 * c4 + i01 * c2 - i03 * c0) * invdet;
    m.i22 = (i30 * s4 - i31 * s2 + i33 * s0) * invdet;
    m.i23 = (-i20 * s4 + i21 * s2 - i23 * s0) * invdet;

    m.i30 = (-i10 * c3 + i11 * c1 - i12 * c0) * invdet;
    m.i31 = (i00 * c3 - i01 * c1 + i02 * c0) * invdet;
    m.i32 = (-i30 * s3 + i31 * s1 - i32 * s0) * invdet;
    m.i33 = (i20 * s3 - i21 * s1 + i22 * s0) * invdet;

    return m;
}

This successfully produced an identity matrix when I multiplied various 3D transformation matrices by the inverse returned from this method. I'm sure you can search/replace to get this into whatever language you'd like.

3
  • 1
    Many thanks for posting, @Robin, this helped me a lot in my C# project. I found one small typo in the code above: in the definition of c5 it should read i31 * i23. After fixing this, matrix inversion works like a charm for me. Mar 8, 2012 at 7:50
  • 1
    Hi @AndersGustafsson, I think you meant the definition of c4 - thanks for the correction - Robin will fix the original.
    – Johnus
    Jun 13, 2013 at 1:10
  • @Johnus You are absolutely right, how silly of me to make this typo when commenting on a typo :-) Thanks for pointing this out. Jun 13, 2013 at 10:34
5

IIRC you can greatly shrink the code and time by precomputing a bunch (12?) 2x2 determinants. Split the matrix in half vertically and compute every 2x2 in both the upper and lower half. One of these smaller determinants is used in every term you'll need for the bigger computation and they each get reused.

Also, don't use a separate determinant function - reuse the sub-determinants you computed for the adjoint to get the determinant.

Oh, just found this.

There are some improvements you can make knowing its a certain kind of transform too.

4
  • Thanks, this saves me a lot of time!
    – Budric
    Apr 12, 2010 at 19:04
  • Very fast, good explanation. Appears to work (haven't run it against a full regression test). Thanks again.
    – Budric
    Apr 12, 2010 at 20:37
  • +1 for the link; however, I think it's a mistake to compute those inverses symbolically... you must realize how many unnecessary multiplications/additions you are performing. It's probably ok as long as this part of the code is not the bottleneck. Apr 13, 2010 at 13:45
  • I use 4x4s for a lot of things, so I prefer the generalized inverse. Like I said, you can do better with specific types of transform. The linked paper is still useful for doing the 3x3 inverse the questioner seems to be using. And you can do even better still if you know the 3x3 is a pure rotation - IIRC it's inverse is the transpose.
    – phkahler
    Apr 14, 2010 at 13:30
2

I believe the only way to compute an inverse is to solve n times the equation: A x = y, where y spans the unit vectors, i.e., the first one is (1,0,0,0), the second is (0,1,0,0), etc.

(Using the cofactors (Cramer's rule) is a bad idea, unless you want a symbolic formula for the inverse.)

Most linear algebra libraries will allow you to solve those linear systems, and even to compute an inverse. Example in python (using numpy):

from numpy.linalg import inv
inv(A) # here you go

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.