# Invert 4x4 Matrix Gaussian Elimination

###
#1
Members - Reputation: **122**

Posted 05 March 2004 - 05:38 AM

###
#4
Members - Reputation: **122**

Posted 05 March 2004 - 09:43 AM

It is hard to mulitply anything by 0 and get 1 for the identity matrix. So this poses a problem for me unless matrices that end up with a 4th row 0 0 0 0 have no determinant.

The problem that I have in coding the algorithm for Gaussian elimination is that everything I find on the internet is that they the pages are mostly theory and not a lot of problem solving. If I could see the Gaussian elimination done on a mock up problem then it would not be that difficult to code. If someone can show me a resource that maybe i missed that works out a sample problem and explains the steps then I could code this algorithm. After I am done then I would be glad to post the source here and also the explanation so that people after me can understand what is going on.

ManBot1: There are NxN matrices classes out there on the internet that would be able to solve your problem. My problem is specific on a 4x4 using Gaussian Elimination. Unless your 9x9 problem will be using the Gaussian Elimination method. If you know a lot about this method or have a resource on it please let me know.

###
#6
Members - Reputation: **1437**

Posted 06 March 2004 - 01:12 AM

Note that transformation matrices (such as rotation matrices)

are orthogonal and thus can be inverted by simply transposing

them.

In most 3D applications this is what you want so you don''t have

to implement general matrix inversion.

###
#7
Anonymous Poster_Anonymous Poster_*
Guests - Reputation:

Posted 06 March 2004 - 01:53 AM

(1/det(A)) * TRANSFORM( A )

I mean 4x4 is pretty small, you''d avoid any branching code (mis-predictions etc). You''d be looking at 56 mults + 23 sub/add instructions to complete the inversion, and I''m pretty sure if you played arround a bit on paper you could reduce the number of mults...

###
#8
Members - Reputation: **122**

Posted 06 March 2004 - 06:44 PM

I was not aware of that for transformation matrices. That may prove to be the quicker route to go but I also want to do this just in case a matrix inversion is necessary on down the road.

Anonymous: I have heard that Gaussian Elimination is a more efficient way of calculating an inverse matrix, although I have not actually been able to test this theory I can not give an exact answer. We will find out when I finish coding this thing.

i will be using an evaluation tool by intel to discover the truth. http://www.gamasutra.com/features/20000131/barad_pfv.htm

at the end of this document is the tool.

On another note if anyone knows the answer to this particular problem on speed, Gaussian Elimination vs. 1/det(A)* Transform (A) Please share with us. It would probably save me some headaches and time.

Thanks.

[edited by - optics on March 7, 2004 7:32:39 AM]

###
#9
Members - Reputation: **1060**

Posted 06 March 2004 - 06:54 PM

(1/det)*adjoint : more calculations, less memory accesses

I can't remember where the turnover point is, but for large matrices the number of calculations required by the (1/det)*adjoint method end up taking more time than the memory accesses for gaussian elimination.

Large matrices : Gaussian elimination is faster

Small matrices : (1/det)*adjoint is faster

I think I remember the (1/det)*adjoint method being faster for 4x4 matrices, but don't quote me on that.

[edited by - joanusdmentia on March 7, 2004 1:55:36 AM]

###
#10
Members - Reputation: **122**

Posted 07 March 2004 - 12:01 AM

http://mathworld.wolfram.com/GaussianElimination.html

What is bugging me is the b

_{1}... b

_{k}part of the equation.

Suppose we have a 4x4 matrix

| e11 e12 e13 e14 |

| e21 e22 e23 e24 |

| e31 e32 e33 e34 |

| e41 e42 e43 e44 |

and we wanted to use Gaussian Elimination to get the matrix in Triangular form

| e11 e12 e13 e14 |

| 0 ~e22 ~e23 ~e24 |

| 0 0 ~e33 ~e34 |

| 0 0 0 ~e44 |

then solve for x

_{1}to x

_{k}. I am having a hard time trying to figure out what b

_{1}, b

_{2}, b

_{3}, b

_{4}is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?

Getting closer to the finish line.

[edited by - optics on March 7, 2004 7:34:06 AM]

[edited by - optics on March 7, 2004 10:55:56 AM]

###
#11
Anonymous Poster_Anonymous Poster_*
Guests - Reputation:

Posted 07 March 2004 - 12:02 AM

quote:Original post by optics

darookie: The matrix is of abitrary size( i.e. 4x4 ) but the values are not.

Not wishing to critisise/be unfair etc., but I''m assuming english isn''t your first language - you seem to be mis-understanding ''arbitrary'' - which means something that could be anything (it hasn''t been specified) of a particular class (for example "an arbitrary distance" would be something that is known to be a distance, but the actual distance isn''t known) - ie an arbitrary sized matrix would be one that could be any size.

I think what you are trying to say is that you are making a matrix class for a matrix with a fixed size of 4 x 4 and the values contained by the matrix are not pre-defined.

[btw, I''m the AP who posted earlier]

You are right that in the general case Gaussian elimination is more efficient, however you are using a very small and fixed sized matrix, so I submit that it is likely that the ''classic'' method of inversion is likely to be faster in particular case. This is mostly to do with the reduced number of memory accesses and the lack of need for loops (with the pivot point etc...). I think the minimum number of mults you can get away with in gaussian elimination for a 4 x 4 matrix would be about 50 anyway, and I''m pretty confident I could eliminate 6 multiplies from the ''classic'' method by eliminating common sub-expressions (you''ll find that several of the 2x2 determinates you calculate on the way to the overall 4x4 determinate are repeated several times through out the operation...calc them once & store them in variables, that will probably be more efficient)

###
#12
Moderators - Reputation: **9793**

Posted 07 March 2004 - 07:33 AM

quote:Original post by optics

then solve for x_{1}to x_{k}. I am having a hard time trying to figure out what b_{1}, b_{2}, b_{3}, b_{4}is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?

Apply the Gaussain to :

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Whe it gets reduced to:

| 1 0 0 0 e''11 e''12 e''13 e''14 |

| 0 1 0 0 e''21 e''22 e''23 e''24 |

| 0 0 1 0 e''31 e''32 e''33 e''34 |

| 0 0 0 1 e''41 e''42 e''43 e''44 |

the e'' 4x4 matrix on the right will be your inverse matrix.

###
#13
Anonymous Poster_Anonymous Poster_*
Guests - Reputation:

Posted 07 March 2004 - 08:58 AM

quote:Original post by optics

then solve for x<sub>1</sub> to x<sub>k</sub>. I am having a hard time trying to figure out what b<sub>1</sub>, b<sub>2</sub>, b<sub>3</sub>, b<sub>4</sub> is suppose to be set to. Do we pull these values out of thin air or is there something I overlooked?quote:

The problem is to solve the system Ax = b

If you calculate the inverse of A, then you can calculate inv(A)b to find x.

You can factor A into upper and lower triangular matrices, to get LUx = b. Then you solve the problem Lc = b, then solve Ux = c. Because the matrices are triangular, the systems can be easily solved with substitution.

You can also attach b as an extra column to A, which is called an augmented matrix, and then use elimination to turn A into an upper triangular matrix, applying the same transformations to b that you do to the rows of A. Then you can solve for x by using the triangular matrix, the new values for b, and substitution.

###
#14
Members - Reputation: **122**

Posted 07 March 2004 - 11:20 AM

1) Could you show me the mathematical solution to obtain I(A

^{-1})?

2)If I post my code would you look over and let me know if I am going in the right direction or am just way off base? I am sure the way that I am doing it is not the most efficient but it seems to be the correct way in my opinion so far.

Then again what do I know?

I thank everyone who has given their help and advice. You guys are really the best.

###
#15
Moderators - Reputation: **9793**

Posted 07 March 2004 - 11:27 AM

quote:Original post by optics

1) Could you show me the mathematical solution to obtain I(A^{-1})?

I''m not sure I understand your question. The identity matrix multiplied by any matrix is the other matrix so I(A

^{-1}) would just be A

^{-1}.

quote:

2)If I post my code would you look over and let me know if I am going in the right direction or am just way off base?

Sure.

###
#16
Members - Reputation: **122**

Posted 07 March 2004 - 11:52 AM

Apply the Gaussain to :

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Whe it gets reduced to:

| 1 0 0 0 e''11 e''12 e''13 e''14 |

| 0 1 0 0 e''21 e''22 e''23 e''24 |

| 0 0 1 0 e''31 e''32 e''33 e''34 |

| 0 0 0 1 e''41 e''42 e''43 e''44 |

2) Here is what I have so far -- be gentle please.

class CMatrix4x4

{

public:

CMatrix4x4(void); // default constructor

CMatrix4x4(double c11, double c12, double c13, double c14,

double c21, double c22, double c23, double c24,

double c31, double c32, double c33, double c34,

double c41, double c42, double c43, double c44); // constructor

~CMatrix4x4(); // destructor

double Determinant(); // determinant of the matrix

CMatrix4x4& Inverse(); // inverse of the matrix

private:

double e11, e12, e13, e14,

e21, e22, e23, e24,

e31, e32, e33, e34,

e41, e42, e43, e44;

};

double CMatrix4x4::Determinant()

{

return e11*e22*e33*e44 - e11*e22*e34*e43 + e11*e23*e34*e42 - e11*e23*e32*e44

+ e11*e24*e32*e43 - e11*e24*e33*e42 - e12*e23*e34*e41 + e12*e23*e31*e44

- e12*e24*e31*e43 + e12*e24*e33*e41 - e12*e21*e33*e44 + e12*e21*e34*e43

+ e13*e24*e31*e42 - e13*e24*e32*e41 + e13*e21*e32*e44 - e13*e21*e34*e42

+ e13*e22*e34*e41 - e13*e22*e31*e44 - e14*e21*e32*e43 + e14*e21*e33*e42

- e14*e22*e33*e41 + e14*e22*e31*e43 - e14*e23*e31*e42 + e14*e23*e32*e41;

}

CMatrix4x4& CMatrix4x4::Inverse()

{

double s1, s2, s3, s4;

double negX;

if(!this->Determinant()) // if determinant is 0.0 then matrix has no inverse

{ return *this; }

else // otherwise there is an inverse and we need to find it

{

/************************** Step 1 of Gaussian Elimination *******************************/

if( e21 == 1 ) // swap row 2 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e21; e12 = e22; e13 = e23; e14 = e24;

e21 = s1; e22 = s2; e23 = s3; e24 = s4;

}

else if( e31 == 1 ) // swap row 3 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e31; e12 = e32; e13 = e33; e14 = e34;

e31 = s1; e32 = s2; e33 = s3; e34 = s4;

}

else if( e41 == 1 ) // swap row 4 with row 1 so that e11 == 1

{ // do Gaussian step 1

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e41; e12 = e42; e13 = e43; e14 = e44;

e41 = s1; e42 = s2; e43 = s3; e44 = s4;

}

// if column 1 does not have a 1 in it then we will do math to make e11 == 1

// sort of hackish but we need to get the code working then we can fix it

//row 2 math

negX = -e21;

e21 = e11/e11*negX + e21;

e22 = e12/e11*negX + e22;

e23 = e13/e11*negX + e23;

e24 = e14/e11*negX + e24;

//row 3 math

negX = -e31;

e31 = e11/e11*negX + e31;

e32 = e12/e11*negX + e32;

e33 = e13/e11*negX + e33;

e34 = e14/e11*negX + e34;

//row 4 math

negX = -e41;

e41 = e11/e11*negX + e41;

e42 = e12/e11*negX + e42;

e43 = e13/e11*negX + e43;

e44 = e14/e11*negX + e44;

/********************** Step 2 of Gaussian Elimination *************************************/

//row 3 math

negX = -e32;

e32 = e22/e22*negX + e32;

e33 = e23/e22*negX + e33;

e34 = e24/e22*negX + e34;

//row 4 math

negX = -e42;

e42 = e22/e22*negX + e42;

e43 = e23/e22*negX + e43;

e44 = e24/e22*negX + e44;

/***************************** Step 3 of Gaussian Elimination *******************************/

//row 4 math

negX = -e43;

e43 = e33/e33*negX + e43;

e44 = e34/e33*negX + e44;

}

return *this;

}

Thanks for the help.

###
#17
Moderators - Reputation: **9793**

Posted 07 March 2004 - 12:18 PM

to here:

| e11 e12 e13 e14 1 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

| 1 0 0 0 e'11 e'12 e'13 e'14 |

| 0 1 0 0 e'21 e'22 e'23 e'24 |

| 0 0 1 0 e'31 e'32 e'33 e'34 |

| 0 0 0 1 e'41 e'42 e'43 e'44 |

We can use this (naive) method. First divide the first row by e11.

| 1 e12/e11 e13/e11 e14/e11 1/e11 0 0 0 |

| e21 e22 e23 e24 0 1 0 0 |

| e31 e32 e33 e34 0 0 1 0 |

| e41 e42 e43 e44 0 0 0 1 |

Then subtract multiples of the first row from each of the second, third and fourth rows.

| 1 e12/e11 e13/e11 e14/e11 1/e11 0 0 0 |

| 0 e*22 e*23 e*24 -e21/e11 1 0 0 |

| 0 e*32 e*33 e*34 -e31/e11 0 1 0 |

| 0 e*42 e*43 e*44 -e31/e11 0 0 1 |

(e*ij) representing calculations I'm too lazy to type out by hand. Now divide the second row by e*22. Then subtract multiples of row two from frows one, three and four. And so on.

For the second part: In your code, since you're bothering to calculate the determinant, there's no need to do Gaussian reduction - just use the standard determinant method of calculating the inverse. Most of the work is already done if you have the determinant.

edit: formatting

[edited by - SiCrane on March 7, 2004 7:18:40 PM]

###
#18
Members - Reputation: **122**

Posted 08 March 2004 - 07:39 AM

Second part: you have to find out of the matrix has a determinant otherwise you will be doing calculations for nothing. If a matrix has a determinant of 0 then there is not an inverse otherwise figure out the inverse. I know this seems backwards to you but I also wanted to code the gaussian elimination method for the knowledge.

I will also end up coding an inverse by way of the determinant method and also by darookie''s suggestion of transposition.

Ok i have done the recode and here is what I have by the math of your suggestion. If its wrong please let me know and I will fix it so that others do not try to copy incorrect code.

class CMatrix4x4

{

public:

CMatrix4x4(void); // default constructor

CMatrix4x4(double c11, double c12, double c13, double c14,

double c21, double c22, double c23, double c24,

double c31, double c32, double c33, double c34,

double c41, double c42, double c43, double c44); // constructor

~CMatrix4x4(); // destructor

Identity(); // identity matrix

double Determinant(); // determinant of the matrix

CMatrix4x4& Inverse(); // inverse of the matrix

private:

double e11, e12, e13, e14,

e21, e22, e23, e24,

e31, e32, e33, e34,

e41, e42, e43, e44;

};

CMatrix4x4::Identity()

{

e11 = e22 = e33 = e44 = 1.0f;

e12 = e13 = e14 = e21 = e23 = e24 =

e31 = e32 = e34 = e41 = e42 = e43 = 0.0f;

}

double CMatrix4x4::Determinant()

{

return e11*e22*e33*e44 - e11*e22*e34*e43 + e11*e23*e34*e42 - e11*e23*e32*e44

+ e11*e24*e32*e43 - e11*e24*e33*e42 - e12*e23*e34*e41 + e12*e23*e31*e44

- e12*e24*e31*e43 + e12*e24*e33*e41 - e12*e21*e33*e44 + e12*e21*e34*e43

+ e13*e24*e31*e42 - e13*e24*e32*e41 + e13*e21*e32*e44 - e13*e21*e34*e42

+ e13*e22*e34*e41 - e13*e22*e31*e44 - e14*e21*e32*e43 + e14*e21*e33*e42

- e14*e22*e33*e41 + e14*e22*e31*e43 - e14*e23*e31*e42 + e14*e23*e32*e41;

}

CMatrix4x4& CMatrix4x4::Inverse()

{

double s1, s2, s3, s4;

double negX, div;

CMatrix4x4 m4x4Inv;

m4x4Inv.Identity();

if(!this->Determinant()) // if determinant is 0.0 then matrix has no inverse

{ return *this; }

else // otherwise there is an inverse and we need to find it

{

/************************** Step 1 of Gaussian Elimination *******************************/

if( e21 == 1 ) // swap row 2 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e21; e12 = e22; e13 = e23; e14 = e24;

e21 = s1; e22 = s2; e23 = s3; e24 = s4;

}

else if( e31 == 1 ) // swap row 3 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e31; e12 = e32; e13 = e33; e14 = e34;

e31 = s1; e32 = s2; e33 = s3; e34 = s4;

}

else if( e41 == 1 ) // swap row 4 with row 1 so that e11 == 1

{

s1 = e11; s2 = e12; s3 = e13; s4 = e14;

e11 = e41; e12 = e42; e13 = e43; e14 = e44;

e41 = s1; e42 = s2; e43 = s3; e44 = s4;

}

//if column 1 does not have a 1 in it then we will do math to make e11 == 1

//sort of hackish but we need to get the code working then we can fix it

// row 1 math

div = e11;

e11 = e11/div;

e12 = e12/div;

e13 = e13/div;

e14 = e14/div;

m4x4Inv.e11 = 1/div;

// row 2 math

negX = -e21;

e21 = e11*negX + e21;

e22 = e12*negX + e22;

e23 = e13*negX + e23;

e24 = e14*negX + e24;

m4x4Inv.e21 = negX/div;

// row 3 math

negX = -e31;

e31 = negX + e31;

e32 = e12*negX + e32;

e33 = e13*negX + e33;

e34 = e14*negX + e34;

m4x4Inv.e31 = negX/div;

// row 4 math

negX = -e41;

e41 = negX + e41;

e42 = e12*negX + e42;

e43 = e13*negX + e43;

e44 = e14*negX + e44;

m4x4Inv.e41 = negX/div;

/********************* Step 2 of Gaussian Elimination *************************************/

div = e22;

// row 2 math

e22 = e22/div;

e23 = e23/div;

e24 = e24/div;

m4x4Inv.e22 = 1/div;

// row 1 math

negX = -e12;

e12 = negX + e12;

e13 = e23*negX + e13;

e14 = e24*negX + e14;

m4x4Inv.e12 = negX/div;

// row 3 math

negX = -e32;

e32 = negX + e32;

e33 = e23*negX + e33;

e34 = e24*negX + e34;

m4x4Inv.e32 = negX/div;

// row 4 math

negX = -e42;

e42 = negX + e42;

e43 = e23*negX + e43;

e44 = e24*negX + e44;

m4x4Inv.e42 = negX/div;

/***************************** Step 3 of Gaussian Elimination *******************************/

// row 3 math

div = e33;

e33 = e33/div;

e34 = e34/div;

m4x4Inv.e33 = 1/div;

//row 1 math

negX = -e13;

e13 = negX + e13;

e14 = e34*negX + e14;

m4x4Inv.e13 = negX/div;

// row 2 math

negX = -e23;

e23 = negX + e23;

e24 = e34*negX + e24;

m4x4Inv.e23 = negX/div;

// row 4 math

negX = -e43;

e43 = negX + e43;

e44 = e34*negX + e44;

m4x4Inv.e43 = negX/div;

/***************************** Step 4 of Gaussian Elimination *******************************/

// row 4 math

div = e44;

e44 = e44/div;

m4x4Inv.e44 = 1/div;

// row 1 math

negX = -e14;

e14 = negX + e14;

m4x4Inv.e14 = negX/div;

// row 2 math

negX = -e24;

e24 = negX + e24;

m4x4Inv.e24 = negX/div;

// row 3 math

negX = -e34;

e34 = negX + e34;

m4x4Inv.e34 = negX/div;

// assign m4x4Inv to the matrix element by element

e11 = m4x4Inv.e11; e12 = m4x4Inv.e12; e13 = m4x4Inv.e13; e14 = m4x4Inv.e14;

e21 = m4x4Inv.e21; e22 = m4x4Inv.e22; e23 = m4x4Inv.e23; e24 = m4x4Inv.e24;

e31 = m4x4Inv.e31; e32 = m4x4Inv.e32; e33 = m4x4Inv.e33; e34 = m4x4Inv.e34;

e41 = m4x4Inv.e41; e42 = m4x4Inv.e42; e43 = m4x4Inv.e43; e44 = m4x4Inv.e44;

}

return *this;

}

###
#19
Moderators - Reputation: **9793**

Posted 08 March 2004 - 09:12 AM

quote:Original post by optics

On the first part, I thank you for that but I also want to ask what you meant by naive method.

For one thing, I didn''t bother to check if e11 was 0 or not before I divided the row by it. Nor did I check to see if e*22 was 0 or not before I divided row 2 by it. Also because of the nature of floating point arithmetic, it''s better, instead of just choosing the number in the first row to divide by, swap the rows so that you divide by the number with the largest value in the first column, etc.

quote:

Second part: you have to find out of the matrix has a determinant otherwise you will be doing calculations for nothing.

Well, since you''re throwing away the determinant calculation, you''re already doing calculations for nothing. If you''re doing Gaussian elimination, then you should just procede with the elimination steps until you get to a point where the you don''t have a valid divisor.

As for the actual code, there are a couple issues with the algorithm implementation.

1: you aren''t checking for divide by zero.

2: when you do the computations for the second, third and fourth rows, you don''t update the previous columns in the inverse half of the augmented matrix.

Also purely, in terms of C++ coding style there are some additional issues.

1: your code would look a lot prettier if you you used std::swap() instead of doing manual swaps.

2: the code would also be a lot cleaner if you create a row vector class to encapsulate some of the matrix math going on in the elimination routine.

3: Some of the function names are misleading. Like the Identity() function. It turns the matrix into an identity matrix, but the name sounds like it might return an identity matrix. Similarly the Inverse() function sounds like it should return the matrix inverse instead of making the current matrix into its inverse.

###
#20
Members - Reputation: **122**

Posted 08 March 2004 - 11:31 AM

Todo for the algorigthm: I know I am not checking for divide by zero and I appreciate you pointing that out to me because knowing me I would forget it in the end. But I also plan on posting the finished code for more critiques.

I also thought it was better if e11 == 1 rather than e11 being the largest value in the column. I have not studied this evaluation at all so I would not know best. If you have a link about this I would be more than willing to read and learn from it. It would also be most appreciated.

Could you show the math for step 2 on updating the previous columns? I have tried figuring this out but it looks to me like I am just multiplying the previous columns by -1.

for example for row 3 in step 2

/** Step 2 of Gaussian Elimination **/

// row 3 math

// div = e22

negX = -e32;

e32 = negX + e32;

e33 = e23*negX + e33;

e34 = e24*negX + e34;

m4x4Inv.e31 = negX/div*m4x4Inv.e31;

m4x4Inv.e32 = negX/div;

should I also be adding m4x4Inv.e

_{ij}in this case m4x4Inv.e31 to the same element? Or am I just way off.

I know it might be asking a lot of you but could you write out the math of that step if only for one update so that I can see it?

Also in the case of a division by zero should I just set those elements to zero?

Thank you for the help on this so far.