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.p...MathMatrix.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, 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_Opt...ung_31125.html

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)

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 THX. It was the missing condition for S.X = 2 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.

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

@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.