Archived

This topic is now archived and is closed to further replies.

BarnyardMessiah

Matrix inversion algorithms

Recommended Posts

How big is the matrix? Do you know anything special about it? If it''s positive semidefinite(and in many cases even if it''s not), you can use the method of conjugate gradients. This works really well if it''s sparse.

If it''s a 4x4, it''s probably faster to just use Cramer''s method.

Otherwise, I guess I''d go with LU or one of those other cool methods i learned in numerical analysis...

Share this post


Link to post
Share on other sites
Just one thing: if you''re trying to get an inverse in order to solve a system of linear equations, it''s just not done in practice. You''re better of whith methods of solving systems. You can also note that getting the inverse equals solving the system AX=I where I is the identity of the same size of A.

Like sjelkjd said, it depends on the size of the matrix and of special characteristics you may know in advance of the matrix. I don''t thing there is an algo that takes less than O(n^3), assuming a "normal matrix" (not sparse and such). I believe there are some strassen variations that take about O(n^2.8) but the matrix has to be kind of big in order to compensate versus a normal Gaussian elimination due to the constants associated. But by using Gaussian elimination, you may need to worry about pivoting. Gaussian elimination is also known by not being stable, wich may be bad for ill-conditioned systems. If numerical precision matters, then i suggest Householder reduction wich is stable.

The O(n^3) can lower if the matrix is sparse or has other characteristics. For this types of matrices there are several methods depending on the characteristics of the matrix. See a linear algebra book for this.

These are all direct methods. There are also iterative methods that take arround O(k*(n^2)), although i''m not into these ones. These are methods like Gauss-Seidel that work with guesses.

And be carefull with crammer rule. Again, like sjelkjd said, if the dimension of the matrix is small it''s a good method. But it''s also an exponential one. Not good for greater dimensions of the matrix.

I hope i have helped.

Share this post


Link to post
Share on other sites
I''m using inverses for inverting the objectspace->worldspace matrix.

I don''t think those are sparse, considering there''s usually a rotation or two involved. I don''t think they''re ill-conditioned either.

It''s mainly a tool for my physics engine, my plan was to use the LU algorithm to invert each row of a 4x4 identity matrix. I''ll give those other ones a look and see...

Share this post


Link to post
Share on other sites
This function work #1 for exactly what you want to do.


CMatrix CMatrix::Invert()
{
CMatrix MatTemp;

if( fabs(this->_44 - 1.0f) > .001f)
this->SetIdentity();
if( fabs(this->_14) > .001f || fabs(this->_24) > .001f || fabs(this->_34) > .001f )
this->SetIdentity();

FLOAT fDetInv = 1.0f / ( this->_11 * ( this->_22 * this->_33 - this->_23 * this->_32 ) -
this->_12 * ( this->_21 * this->_33 - this->_23 * this->_31 ) +
this->_13 * ( this->_21 * this->_32 - this->_22 * this->_31 ) );

MatTemp._11 = fDetInv * ( this->_22 * this->_33 - this->_23 * this->_32 );
MatTemp._12 = -fDetInv * ( this->_12 * this->_33 - this->_13 * this->_32 );
MatTemp._13 = fDetInv * ( this->_12 * this->_23 - this->_13 * this->_22 );
MatTemp._14 = 0.0f;

MatTemp._21 = -fDetInv * ( this->_21 * this->_33 - this->_23 * this->_31 );
MatTemp._22 = fDetInv * ( this->_11 * this->_33 - this->_13 * this->_31 );
MatTemp._23 = -fDetInv * ( this->_11 * this->_23 - this->_13 * this->_21 );
MatTemp._24 = 0.0f;

MatTemp._31 = fDetInv * ( this->_21 * this->_32 - this->_22 * this->_31 );
MatTemp._32 = -fDetInv * ( this->_11 * this->_32 - this->_12 * this->_31 );
MatTemp._33 = fDetInv * ( this->_11 * this->_22 - this->_12 * this->_21 );
MatTemp._34 = 0.0f;

MatTemp._41 = -( this->_41 * MatTemp._11 + this->_42 * MatTemp._21 + this->_43 * MatTemp._31 );
MatTemp._42 = -( this->_41 * MatTemp._12 + this->_42 * MatTemp._22 + this->_43 * MatTemp._32 );
MatTemp._43 = -( this->_41 * MatTemp._13 + this->_42 * MatTemp._23 + this->_43 * MatTemp._33 );
MatTemp._44 = 1.0f;

return MatTemp;
}


Fantasio

[edited by - fantasio on November 16, 2003 2:29:43 AM]

Share this post


Link to post
Share on other sites
v71 is right. If a matrix is orthogonal (all its vectors 90 degrees to each other) and normal (all its vectors are unit vectors), it''s transpose is equivalent to its inverse.

With objectspace->worldspace matrix conversions this is generally the case, unless you are doing wierd scaling of your object into worldspace.

Of course, Fantasio''s method will probably work for any invertible matrix (it does not do any invertible-tests and will cause divide by zero in non-invertible cases).

Share this post


Link to post
Share on other sites
Yeah, I implemented 1/det(A)*adjoint(A) yesterday (basically what Fantasio posted), but it didn''t occure to me that just transposing an orthogonal unit matrix was the same as inverting it.

Makes sence since you''ve explained it.

But I''ve also implemented a transpose function in my Matrix class so I guess I have the best of both worlds.

Share this post


Link to post
Share on other sites