# Fastest way to inverse an orthogonal 4x4 matrix?

This topic is 5408 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

I know about the Kramer's rule (and have that working), the Gaussian method, and know that the inverse of a rotational 3x3 matrix is it's transpose. But I think I'm missing something. Since the matrix is orthogonal and I can calculate the 3x3 part of it by transposing it, shouldn't there be some easier method to calculate the remaining 3 numbers (the translate part)? I heard somewhere that they should be flipped from positive to negative, but that's not quite working. What am I missing? Here's an example... I have this matrix:
0  4  8  12
1  5  9  13
2  6  10 14
3  7  11 15

X:	0.819300	-0.196000	0.538800	-0.591540
Y:	0.301200	0.946800	-0.113700	-1.476260
Z:	-0.487900	0.255500	0.834700	-0.163410
W:	0.000000	0.000000	0.000000	1.000000

Using the Kramer's rule, the inverse of it is:
X:	0.819302	0.301249	-0.487825	0.849656
Y:	-0.195927	0.946701	0.255428	1.323417
Z:	0.538872	-0.113697	0.834705	0.287318
W:	0.000000	0.000000	0.000000	1.000000

But the center can be found by simply transposing the 3x3 part:
X:	0.819300	0.301200	-0.487900	-0.591540
Y:	-0.196000	0.946800	0.255500	-1.476260
Z:	0.538800	-0.113700	0.834700	-0.163410
W:	0.000000	0.000000	0.000000	1.000000

So the question is, how to calculate the m[12], m[13] and m[14] quickly? I suppose simply taking a part of the code from Kramer's rule for those 3 would do it, but that's cheap... Does anyone know a better way?

##### Share on other sites
if you have 4x4 matrix A that acts as rotation matrix3x3 B with some translation V ,
that is, forward transform for X to X' is X'=A*X=B*X+V

So backward transform will be (matrix inverse is denoted by ~)
X = (~A)*X' = B'*X'+V' = (~B)*(X'-V) = (~B)*X'-(~B)*V
so
B' = (~B)
V' = - B'*V

that is, to compute inverse of 4X4 matrix with rotation and translation only(translation part is a last column) you need to:

transpose 3x3 matrix in top-left corner.
Multiply translation part by 3x3 matrix in top-left corner.
Negate translation part.

##### Share on other sites
also, note that the lower row has three zeros and a one in it. if you know how to invert matrices, that should be enough to simplify things.

##### Share on other sites
"3x3 part of it by transposing it"
/!\ Orthonormal (no scaling) is not orthogonal.

That's obvious if you use the linear properties.
Decompose your matrix by blocks as Mt*Mr. You first rotate by the 3x3 submatrix R then translate by T. Thus the inverse is
M^-1= Mr^-1 * Mt^-1
Consider R is made of columns X, Y, Z (the local axes).

 Mt             Mr            M 1  0  0  Tx    .  .  .  0    Xx Yx Zx Tx 0  1  0  Ty    .  R  .  0    Xy Yy Zy Ty 0  0  1  Tz    .  .  .  0    Xz Yz Zz Tz 0  0  0  1     0  0  0  1    0  0  0  1 Mt^-1          Mr^-1        M^-1 1  0  0 -Tx    .  .  .  0   Xx Xy Xz -X*Tx 0  1  0 -Ty    . R^-1.  0   Yx Yy Yz -Y*Ty 0  0  1 -Tz    .  .  .  0   Zx Zy Zz -Z*Tz 0  0  0  1     0  0  0  1   0  0  0  1

Comment :
In most typical cases, I would not even dare to inverse such a matrix. I would write a special routine named "transform vector by transposed" :

P' <- { X*(P-T), Y*(P-T), Z*(P-T) }

Edit : Dmytry and Eelco. Always the same club. You both beat me by a few seconds :) LOL !!!

##### Share on other sites
Thanks! This works:

// Inverses an orthogonal matrixvoid CMatrix4x4::InvertOrthogonal(){	float fOut[16];	memcpy(fOut, fMat, sizeof(float) * 16);	fMat[1]  =  fOut[4];	fMat[2]  =  fOut[8];	fMat[4]  =  fOut[1];	fMat[6]  =  fOut[9];	fMat[8]  =  fOut[2];	fMat[9]  =  fOut[6];	fMat[12] = -(fOut[12] * fMat[0]  +  fOut[13] * fMat[4]  +  fOut[14] * fMat[8]);	fMat[13] = -(fOut[12] * fMat[1]  +  fOut[13] * fMat[5]  +  fOut[14] * fMat[9]);	fMat[14] = -(fOut[12] * fMat[2]  +  fOut[13] * fMat[6]  +  fOut[14] * fMat[10]);}

##### Share on other sites
Quote:
 Original post by Charles B"3x3 part of it by transposing it"/!\ Orthonormal (no scaling) is not orthogonal.

Hehe. I'm a math noob, I can make such mistakes. x_X

At least it starts with the same "Ortho"! :D
Quote:
 That's obvious if you use the linear properties.Decompose your matrix by blocks as Mt*Mr. You first rotate by the 3x3 submatrix R then translate by T. Thus the inverse isM^-1= Mr^-1 * Mt^-1Consider R is made of columns X, Y, Z (the local axes). Mt Mr M 1 0 0 Tx . . . 0 Xx Yx Zx Tx 0 1 0 Ty . R . 0 Xy Yy Zy Ty 0 0 1 Tz . . . 0 Xz Yz Zz Tz 0 0 0 1 0 0 0 1 0 0 0 1 Mt^-1 Mr^-1 M^-1 1 0 0 -Tx . . . 0 Xx Xy Xz -X*Tx 0 1 0 -Ty . R^-1. 0 Yx Yy Yz -Y*Ty 0 0 1 -Tz . . . 0 Zx Zy Zz -Z*Tz 0 0 0 1 0 0 0 1 0 0 0 1Comment :In most typical cases, I would not even dare to inverse such a matrix. I would write a special routine named "transform vector by transposed" :P' <- { X*(P-T), Y*(P-T), Z*(P-T) }Edit : Dmytry and Eelco. Always the same club. You both beat me by a few seconds :) LOL !!!

You guys are all awesome. ;) Much appreciated~

##### Share on other sites
Quote:
 Original post by MichaelMookThanks! This works:

OK but pls, optimize it a bit now that it works.

Use temporary floats instead of copying the full matrix on the stack. Also beware memcpy is sometimes very slow, in debug for instance. Even John Carmack fell in the trap once. He wondered why his code was so slow and it took him several hours to see it came from the std C library implementation. It's somthing I knew before he wrote his first books, I have always checked what the compilers do. I am a speedy code addict, low level biased :)

You need 9 temporary variables. With some luck, since the FPU has 8 registers, only one such variable will be on the stack.

Write something like :

float f4=M[4], f8=M[8], etc ...;

M[1]=f4;

##### Share on other sites
Obviously for transformation matrices the other presentations here, especially Charles B's, are THE appropriate way to do this if you absolutely must have the inverse.

But, suppose you have a general, non-transformation matrix...

Just a note. For matrices around 4x4 and above, Cramer's rule is going to be the slowest method! It's an O((n+1)!)) operation. For a general matrix 3x3 or larger, use Gaussian Elimination with partial pivoting or one of the more advanced direct method. Very large matrices (maybe 1000x1000 or larger) require an iterative method, since roundoff will kill Gaussian Elimination.

[edited by grhodes_at_work to correct the operation count for Cramer's Rule. Its far, far worse than O(n3).]

[Edited by - grhodes_at_work on September 25, 2004 7:32:48 PM]

##### Share on other sites
Quote:
Original post by Charles B
Quote:
 Original post by MichaelMookThanks! This works:

OK but pls, optimize it a bit now that it works.

Use temporary floats instead of copying the full matrix on the stack. Also beware memcpy is sometimes very slow, in debug for instance. Even John Carmack fell in the trap once. He wondered why his code was so slow and it took him several hours to see it came from the std C library implementation. It's somthing I knew before he wrote his first books, I have always checked what the compilers do. I am a speedy code addict, low level biased :)

You need 9 temporary variables. With some luck, since the FPU has 8 registers, only one such variable will be on the stack.

Write something like :

float f4=M[4], f8=M[8], etc ...;

M[1]=f4;

Thanks, I'll keep that in mind. I didn't notice any change in framerate, but I will keep it in mind for the future when I have more happening in the scene and will be able to judge more precisely.

Quote:
 Just a note. For matrices around 4x4 and above, Cramer's rule is going to be the slowest method!

Yup, I saw your other posts ;). I'm not really using it though. I also ended up making two other versions of the above function, as it seems I was getting eventual rounding errors by not using the determinant. And the other one merely flips the translation coordinates (seems to work nicely for billboarding of the sun and requires the least number of calculations).

##### Share on other sites
Sorry i haven't read the whole thread but i've been using Cramer's rule expressed differently to get the inverse matrix explicitly for 4x4 matrix and below:

M-1 = (1 / det(M)) * adj(M)

where det = determinate &

i thought this was the fastest way for matrices upto 4x4? anything bigger than 4x4 use Gaussian elimination.

and for orthogonal matrices:

M-1 = MT

[Edited by - snk_kid on September 24, 2004 5:40:23 PM]

• 21
• 14
• 9
• 17
• 13