View Full Version : Math: Calculation of inverse Matrizes
BenBE
14-10-2004, 11:20 AM
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
14-10-2004, 01:46 PM
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<sup>-1</sup> = (1/|A|) * (adjoint of A)
(If |A| = 0, then there is no inverse matrix.)
BenBE
14-10-2004, 06:52 PM
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
14-10-2004, 07:17 PM
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
20-10-2004, 02:00 PM
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
20-10-2004, 06:18 PM
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
21-10-2004, 12:50 PM
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
21-10-2004, 09:08 PM
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
22-10-2004, 01:41 PM
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[0][0] <> 0 then
Result[0] := ARightSide[0] / ALeftSide[0][0] else
ShowMessage('Equation is unsolvable - it is inconsistent');
Exit;
end;
// solve the denominator of Cramer's equation
// only needs to be calculated once
LDenominator := Determinant(ALeftSide);
if LDenominator = 0 then begin
ShowMessage('Equations are unsolvable - they are either dependant or inconsistent');
Exit;
end;
// calculate the inverse so that the divide is only done once
LDenominatorInv := 1 / LDenominator;
// create numerator for first variable
SetLength(LNumeratorDet, LOrder);
for LRow := 0 to LOrder - 1 do begin
SetLength(LNumeratorDet[LRow], LOrder);
LNumeratorDet[LRow][0] := ARightSide[LRow];
for LCol := 1 to LOrder - 1 do begin
LNumeratorDet[LRow][LCol] := ALeftSide[LRow][LCol];
end;
end;
// solve the first variable
Result[0] := Determinant(LNumeratorDet) * LDenominatorInv;
// solve the remaining variables
for LCol := 1 to LOrder - 1 do begin
// restore the previous column with the left side values
for LRow := 0 to LOrder - 1 do begin
LNumeratorDet[LRow][LCol-1] := ALeftSide[LRow][LCol-1];
end;
// replace the current column with the right side values
for LRow := 0 to LOrder - 1 do begin
LNumeratorDet[LRow][LCol] := ARightSide[LRow];
end;
// solve for this variable
Result[LCol] := Determinant(LNumeratorDet) * LDenominatorInv;
end;
end;
BenBE
22-10-2004, 05:18 PM
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
25-10-2004, 09:07 PM
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?
Powered by vBulletin® Version 4.2.5 Copyright © 2024 vBulletin Solutions Inc. All rights reserved.