Math: Calculation of inverse Matrizes

[BenBE]

I've been wondering if there's a way to make calculating inverse-Matrizzes faster then by usings the Gauss Algorith and apply it's opperations onto a identity matrix.

If anyone has an idea or could give me some information on other algorithms available for this task I'd appreciate to hear about them.

(On request I can send my current solution I don't post it now as it is dependen of a lot of code)

FYI: It's function MatrixInverse in http://cvs.sourceforge.net/viewcvs.py/omorphia/omorphia/maths/OMathMatrix.pas

[Useless Hacker]

This is the only method I know of for finding an inverse matrix.

To find the inverse of a matrix A:
1) Find the determinant of A, |A|.

2) Find the adjoint of A, which is the transpose of the matrix of cofactors of A.

3) The inverse of A is given by:
A^-1 = (1/|A|) * (adjoint of A)

(If |A| = 0, then there is no inverse matrix.)

[BenBE]

Sorry for the dumb question, but I'm new to matrix calculations:
What are the cofactors of a matrix? And how to get them?

The limitation at the bottom is clear if you look at the definition ;-)

THX so far.

[Useless Hacker]

A quick search on Google comes up with this PDF file (http://www.mathcentre.ac.uk/resources/leaflets/firstaidkits/5_5.pdf), which explains the whole process better than I can ;).

[peterbone]

here's a function I wrote to calculate a determinant of any order recursively. Just in case you needed one.


TFloatArray = array of Extended;
T2DFloatArray = array of TFloatArray;

// solve a determinant of any order
function Determinant(ADet : T2DFloatArray) : Extended;
Var
LCol, LRow, LOrder, LNewOrder : Integer;
LSign : boolean;
LMinor : T2DFloatArray;
begin
LOrder := length(ADet);

if LOrder = 2 then begin

// return the result
Result := ADet[0][0] * ADet[1][1] - ADet[0][1] * ADet[1][0];

end else begin

// reduce the determinant by minors and solve determinants recursively
LNewOrder := LOrder - 1;

// create first minor
SetLength(LMinor, LNewOrder);
for LRow := 0 to LNewOrder - 1 do begin
SetLength(LMinor[LRow], LNewOrder);
for LCol := 0 to LNewOrder - 1 do
LMinor[LRow][LCol] := ADet[LRow+1][LCol+1];
end;

// solve first minor
Result := ADet[0][0] * Determinant(LMinor);

LSign := False; // negative

// solve remaining minors and accumulate result
for LCol := 1 to LNewOrder do begin

// create the minor for this column by modifying the previous minor
for LRow := 1 to LNewOrder do
LMinor[LRow-1][LCol-1] := ADet[LRow][LCol-1];

// accumulate terms to get result
if LSign then Result := Result + ADet[0][LCol] * Determinant(LMinor)
else Result := Result - ADet[0][LCol] * Determinant(LMinor);

LSign := not LSign; // alternate sign of terms

end; // minor columns loop

end;
end;

[BenBE]

For Determinants I already wrote a very fast implementation (0.95 ms per 1000x1000 Matrix) which you can find here: http://www.delphi-forum.de/topic_Optimierung+DeterminantenBerechnung_31125.ht ml

For Inverse Matrizes I had no time yet to implement the suggestions given in the PDF file which helped a lot.

I'll tell when I've something that works (or not) :D

THX NE-Way

[peterbone]

ok, your determinant function looks much faster than mine. What method does it use? I tried to use it but I couldn't get it to work - I copied it exactly as it is in the link you gave. For a 2x2 matrix with values ((1 2)(3 4)) it should give -2 but gives 0. What could I be doing wrong?

BTW, I think there's a numerical solution to find the inverse of a matrix, which is fastest for large matrices.

Thanks

Peter Bone

[BenBE]

It uses the mathematical definition of finding the determinant, i.e. it calculates


a1 b1 c1
a2 b2 c2
a3 b3 c3


Det =
//diagonals from left top to right bottom
a1 * b2 * c3 +
b1 * c2 * a3 +
c1 * a2 * b3
-
(
//diagonals from right top to left bottom
c1 * b2 * a3 +
b1 * a2 * c3 +
a1 * c2 * b3
);

For 2x2 Matrizes I already noticed this problem. It's a problem with the mathematical solution I used. For 2x2 Matrizes it simplifies to Det := 0, which is wrong.

But your routine gave me the solution for the problem :D THX. It was the missing condition for S.X = 2 :D I'll include it into my source and it should be right.

But the main problem remains ... I'll tell, when I've something that works. But could take some days.

[peterbone]

ok, could you post your determinant code when you've fixed it please.
Why do you need to calculate a matrix inverse? If you're solving a system of linear equations then you don't have to calculate the inverse - you can work it out quicker using Cramers method. Here's my code:

// solve a system of linear equations of any order using Cramer's method
function SolveLinearEquations(ALeftSide : T2DFloatArray ; ARightSide : TFloatArray) : TFloatArray;
Var
LDenominator, LDenominatorInv : Extended;
LOrder : integer;
LRow, LCol : integer;
LNumeratorDet : T2DFloatArray;
begin
LOrder := length(ARightSide);
SetLength(Result, LOrder);

// special case - order = 1
if LOrder = 1 then begin
if ALeftSide&#91;0&#93;&#91;0&#93; <> 0 then
Result&#91;0&#93; &#58;= ARightSide&#91;0&#93; / ALeftSide&#91;0&#93;&#91;0&#93; else
ShowMessage&#40;'Equation is unsolvable - it is inconsistent'&#41;;
Exit;
end;

// solve the denominator of Cramer's equation
// only needs to be calculated once
LDenominator &#58;= Determinant&#40;ALeftSide&#41;;

if LDenominator = 0 then begin
ShowMessage&#40;'Equations are unsolvable - they are either dependant or inconsistent'&#41;;
Exit;
end;

// calculate the inverse so that the divide is only done once
LDenominatorInv &#58;= 1 / LDenominator;

// create numerator for first variable
SetLength&#40;LNumeratorDet, LOrder&#41;;
for LRow &#58;= 0 to LOrder - 1 do begin
SetLength&#40;LNumeratorDet&#91;LRow&#93;, LOrder&#41;;
LNumeratorDet&#91;LRow&#93;&#91;0&#93; &#58;= ARightSide&#91;LRow&#93;;
for LCol &#58;= 1 to LOrder - 1 do begin
LNumeratorDet&#91;LRow&#93;&#91;LCol&#93; &#58;= ALeftSide&#91;LRow&#93;&#91;LCol&#93;;
end;
end;

// solve the first variable
Result&#91;0&#93; &#58;= Determinant&#40;LNumeratorDet&#41; * LDenominatorInv;

// solve the remaining variables
for LCol &#58;= 1 to LOrder - 1 do begin

// restore the previous column with the left side values
for LRow &#58;= 0 to LOrder - 1 do begin
LNumeratorDet&#91;LRow&#93;&#91;LCol-1&#93; &#58;= ALeftSide&#91;LRow&#93;&#91;LCol-1&#93;;
end;

// replace the current column with the right side values
for LRow &#58;= 0 to LOrder - 1 do begin
LNumeratorDet&#91;LRow&#93;&#91;LCol&#93; &#58;= ARightSide&#91;LRow&#93;;
end;

// solve for this variable
Result&#91;LCol&#93; &#58;= Determinant&#40;LNumeratorDet&#41; * LDenominatorInv;
end;
end;

[BenBE]

I need the inverse of an Matrix for a maths library I'm writing related to the Project Omorphia (http://sf.net/project/omorphia/).

I know that solving a linear equation system doesn't require to know it's inverse matrix, but the steps of transforming an matrix A into an identity matrix can be used, when applied to an identitiy matrix to create the inverse of the matrix A. But this algorithm seems to be rather slow and a nice place for various bugs :D

@My MatrixDeterminant-Routine: I tested it today with the Delphi Random-Function and found that it's not very accurate (or my tested MatrixInverse routine wasn't). The updated source can be found under the previously given link.

[BenBE]

This is my currently used function for calculating the inverse Matrix. I hope it's complete. For the MatrixDeterminant-Routine see the link I gave in the first post.

Function MatrixInverse(Const M: TMatrix): TMatrix;
Var
S: TMatrixSize;
// X, Y: Integer;
// I, J, K: Integer;
// Swapped: Boolean;
// Temp: TMatrix;
Det: Extended;
Begin
//Get the dimensions of the matrix
S := GetMatrixDim(M);

//Check if both Matrizes have at least one col and row
If (S.X <= 0) Or (S.Y <= 0) Then
Raise EMathError.Create(FmtInvalidDim);
If S.X <> S.Y Then
Raise EMathError.Create(FmtInvalidDim);

Det := MatrixDeterminant(M);
If Det = 0 Then
Raise EMathError.Create(FmtMatrixNotInversable);

Result := MatrixRMul(MatrixAdjoint(M), 1 / Det);
End;

Function MatrixAdjoint(Const M: TMatrix): TMatrix;
Var
S: TMatrixSize;
T: TMatrix;

Function MinorAndCofactor(X, Y: Integer): Extended;
Var
I, J: Integer;
A, B: Integer;
C: TMatrix;
Begin
C := SetupMatrix(S.Y - 1, S.X - 1);
B := 0;
For J := 0 To S.Y - 2 Do
Begin
If B = Y Then
Inc(B);

A := 0;
For I := 0 To S.X - 2 Do
Begin
If A = X Then
Inc(A);
C[J, I] := T[B, A];
Inc(A);
End;

Inc(B);
End;

Asm
// Result := MatrixDeterminant(C);
MOV EAX, DWORD PTR [C]
CALL MatrixDeterminant
//Apply the cofactor sign change rule
MOV EAX, DWORD PTR [X]
XOR EAX, DWORD PTR [Y]
AND EAX, $00000001
JZ @@Finish
//Change the sign
FCHS
@@Finish:
FSTP Result
End;
End;

Var
X, Y: Integer;
Begin
S := GetMatrixDim(M);

If (S.X = 0) Or (S.X <> S.Y) Then
Raise EMathError.Create(FmtInvalidDim);

T := MatrixTranspose(M);

Result := SetupMatrix(S);

For Y := 0 To S.Y - 1 Do
For X := 0 To S.X - 1 Do
Result[Y, X] := MinorAndCofactor(X, Y);
End;

Any improvement possible?