Simple 3x3 matrix inverse code (C++)

Question

What's the easiest way to compute a 3x3 matrix inverse?

I'm just looking for a short code snippet that'll do the trick for non-singular matrices, possibly using Cramer's rule. It doesn't need to be highly optimized. I'd prefer simplicity over speed. I'd rather not link in additional libraries.

Answer 1

With all due respect to our unknown (yahoo) poster, I look at code like that and just die a little inside. Alphabet soup is just so insanely difficult to debug. A single typo anywhere in there can really ruin your whole day. Sadly, this particular example lacked variables with underscores. It's so much more fun when we have a_b-c_d*e_f-g_h. Especially when using a font where _ and - have the same pixel length.

Taking up Suvesh Pratapa on his suggestion, I note:

Given 3x3 matrix:
       y0x0  y0x1  y0x2
       y1x0  y1x1  y1x2
       y2x0  y2x1  y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];

(A) When taking a minor of a 3x3 array, we have 4 values of interest. The lower X/Y index is always 0 or 1. The higher X/Y index is always 1 or 2. Always! Therefore:

double determinantOfMinor( int          theRowHeightY,
                           int          theColumnWidthX,
                           const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  int x1 = theColumnWidthX == 0 ? 1 : 0;  /* always either 0 or 1 */
  int x2 = theColumnWidthX == 2 ? 1 : 2;  /* always either 1 or 2 */
  int y1 = theRowHeightY   == 0 ? 1 : 0;  /* always either 0 or 1 */
  int y2 = theRowHeightY   == 2 ? 1 : 2;  /* always either 1 or 2 */

  return ( theMatrix [y1] [x1]  *  theMatrix [y2] [x2] )
      -  ( theMatrix [y1] [x2]  *  theMatrix [y2] [x1] );
}

(B) Determinant is now: (Note the minus sign!)

double determinant( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  return ( theMatrix [0] [0]  *  determinantOfMinor( 0, 0, theMatrix ) )
      -  ( theMatrix [0] [1]  *  determinantOfMinor( 0, 1, theMatrix ) )
      +  ( theMatrix [0] [2]  *  determinantOfMinor( 0, 2, theMatrix ) );
}

(C) And the inverse is now:

bool inverse( const double theMatrix [/*Y=*/3] [/*X=*/3],
                    double theOutput [/*Y=*/3] [/*X=*/3] )
{
  double det = determinant( theMatrix );

    /* Arbitrary for now.  This should be something nicer... */
  if ( ABS(det) < 1e-2 )
  {
    memset( theOutput, 0, sizeof theOutput );
    return false;
  }

  double oneOverDeterminant = 1.0 / det;

  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
        /* Rule is inverse = 1/det * minor of the TRANSPOSE matrix.  *
         * Note (y,x) becomes (x,y) INTENTIONALLY here!              */
      theOutput [y] [x]
        = determinantOfMinor( x, y, theMatrix ) * oneOverDeterminant;

        /* (y0,x1)  (y1,x0)  (y1,x2)  and (y2,x1)  all need to be negated. */
      if( 1 == ((x + y) % 2) )
        theOutput [y] [x] = - theOutput [y] [x];
    }

  return true;
}

And round it out with a little lower-quality testing code:

void printMatrix( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  for ( int y = 0;  y < 3;  y ++ )
  {
    cout << "[  ";
    for ( int x = 0;  x < 3;  x ++   )
      cout << theMatrix [y] [x] << "  ";
    cout << "]" << endl;
  }
  cout << endl;
}

void matrixMultiply(  const double theMatrixA [/*Y=*/3] [/*X=*/3],
                      const double theMatrixB [/*Y=*/3] [/*X=*/3],
                            double theOutput  [/*Y=*/3] [/*X=*/3]  )
{
  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
      theOutput [y] [x] = 0;
      for ( int i = 0;  i < 3;  i ++ )
        theOutput [y] [x] +=  theMatrixA [y] [i] * theMatrixB [i] [x];
    }
}

int
main(int argc, char **argv)
{
  if ( argc > 1 )
    SRANDOM( atoi( argv[1] ) );

  double m[3][3] = { { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) } };
  double o[3][3], mm[3][3];

  if ( argc <= 2 )
    cout << fixed << setprecision(3);

  printMatrix(m);
  cout << endl << endl;

  SHOW( determinant(m) );
  cout << endl << endl;

  BOUT( inverse(m, o) );
  printMatrix(m);
  printMatrix(o);
  cout << endl << endl;

  matrixMultiply (m, o, mm );
  printMatrix(m);
  printMatrix(o);
  printMatrix(mm);  
  cout << endl << endl;
}

---

Afterthought:

You may also want to detect very large determinants as round-off errors will affect your accuracy!

Answer 2

I went ahead and wrote it in python since I think it's a lot more readable than in c++ for a problem like this. The function order is in order of operations for solving this by hand via this video. Just import this and call "print_invert" on your matrix.

def print_invert (matrix):
  i_matrix = invert_matrix (matrix)
  for line in i_matrix:
    print (line)
  return

def invert_matrix (matrix):
  determinant = str (determinant_of_3x3 (matrix))
  cofactor = make_cofactor (matrix)
  trans_matrix = transpose_cofactor (cofactor)

  trans_matrix[:] = [[str (element) +'/'+ determinant for element in row] for row in trans_matrix]

  return trans_matrix

def determinant_of_3x3 (matrix):
  multiplication = 1
  neg_multiplication = 1
  total = 0
  for start_column in range (3):
    for row in range (3):
      multiplication *= matrix[row][(start_column+row)%3]
      neg_multiplication *= matrix[row][(start_column-row)%3]
    total += multiplication - neg_multiplication
    multiplication = neg_multiplication = 1
  if total == 0:
    total = 1
  return total

def make_cofactor (matrix):
  cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  matrix_2x2 = [[0,0],[0,0]]
  # For each element in matrix...
  for row in range (3):
    for column in range (3):

      # ...make the 2x2 matrix in this inner loop
      matrix_2x2 = make_2x2_from_spot_in_3x3 (row, column, matrix)
      cofactor[row][column] = determinant_of_2x2 (matrix_2x2)

  return flip_signs (cofactor)

def make_2x2_from_spot_in_3x3 (row, column, matrix):
  c_count = 0
  r_count = 0
  matrix_2x2 = [[0,0],[0,0]]
  # ...make the 2x2 matrix in this inner loop
  for inner_row in range (3):
    for inner_column in range (3):
      if row is not inner_row and inner_column is not column:
        matrix_2x2[r_count % 2][c_count % 2] = matrix[inner_row][inner_column]
        c_count += 1
    if row is not inner_row:
      r_count += 1
  return matrix_2x2

def determinant_of_2x2 (matrix):
  total = matrix[0][0] * matrix [1][1]
  return total - (matrix [1][0] * matrix [0][1])

def flip_signs (cofactor):
  sign_pos = True 
  # For each element in matrix...
  for row in range (3):
    for column in range (3):
      if sign_pos:
        sign_pos = False
      else:
        cofactor[row][column] *= -1
        sign_pos = True
  return cofactor

def transpose_cofactor (cofactor):
  new_cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  for row in range (3):
    for column in range (3):
      new_cofactor[column][row] = cofactor[row][column]
  return new_cofactor

Answer 3

//Function for inverse of the input square matrix 'J' of dimension 'dim':

vector<vector<double > > inverseVec33(vector<vector<double > > J, int dim)
{
//Matrix of Minors
 vector<vector<double > > invJ(dim,vector<double > (dim));
for(int i=0; i<dim; i++)
{
    for(int j=0; j<dim; j++)
    {
        invJ[i][j] = (J[(i+1)%dim][(j+1)%dim]*J[(i+2)%dim][(j+2)%dim] -
                      J[(i+2)%dim][(j+1)%dim]*J[(i+1)%dim][(j+2)%dim]);
    }
}

//determinant of the matrix:
double detJ = 0.0;
for(int j=0; j<dim; j++)
{ detJ += J[0][j]*invJ[0][j];}

//Inverse of the given matrix.
 vector<vector<double > > invJT(dim,vector<double > (dim));
 for(int i=0; i<dim; i++)
{
    for(int j=0; j<dim; j++)
    {
        invJT[i][j] = invJ[j][i]/detJ;
    }
}

return invJT;
}

void main()
{
    //given matrix:
vector<vector<double > > Jac(3,vector<double > (3));
Jac[0][0] = 1; Jac[0][1] = 2;  Jac[0][2] = 6;
Jac[1][0] = -3; Jac[1][1] = 4;  Jac[1][2] = 3;
Jac[2][0] = 5; Jac[2][1] = 1;  Jac[2][2] = -4;`

//Inverse of the matrix Jac:
vector<vector<double > > JacI(3,vector<double > (3));
    //call function and store inverse of J as JacI:
JacI = inverseVec33(Jac,3);
}

Answer 4

//Title: Matrix Header File
//Writer: Say OL
//This is a beginner code not an expert one
//No responsibilty for any errors
//Use for your own risk
using namespace std;
int row,col,Row,Col;
double Coefficient;
//Input Matrix
void Input(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
        {
            cout<<"e["<<row<<"]["<<col<<"]=";
            cin>>Matrix[row][col];
        }
}
//Output Matrix
void Output(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
    {
        for(col=1;col<=Col;col++)
            cout<<Matrix[row][col]<<"\t";
        cout<<endl;
    }
}
//Copy Pointer to Matrix
void CopyPointer(double (*Pointer)[9],double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            Matrix[row][col]=Pointer[row][col];
}
//Copy Matrix to Matrix
void CopyMatrix(double MatrixInput[9][9],double MatrixTarget[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTarget[row][col]=MatrixInput[row][col];
}
//Transpose of Matrix
double MatrixTran[9][9];
double (*(Transpose)(double MatrixInput[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTran[col][row]=MatrixInput[row][col];
    return MatrixTran;
}
//Matrix Addition
double MatrixAdd[9][9];
double (*(Addition)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixAdd[row][col]=MatrixA[row][col]+MatrixB[row][col];
    return MatrixAdd;
}
//Matrix Subtraction
double MatrixSub[9][9];
double (*(Subtraction)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixSub[row][col]=MatrixA[row][col]-MatrixB[row][col];
    return MatrixSub;
}
//Matrix Multiplication
int mRow,nCol,pCol,kcol;
double MatrixMult[9][9];
double (*(Multiplication)(double MatrixA[9][9],double MatrixB[9][9],int mRow,int nCol,int pCol))[9]
{
    for(row=1;row<=mRow;row++)
        for(col=1;col<=pCol;col++)
        {
            MatrixMult[row][col]=0.0;
            for(kcol=1;kcol<=nCol;kcol++)
                MatrixMult[row][col]+=MatrixA[row][kcol]*MatrixB[kcol][col];
        }
    return MatrixMult;
}
//Interchange Two Rows
double RowTemp[9][9];
double MatrixInter[9][9];
double (*(InterchangeRow)(double MatrixInput[9][9],int Row,int Col,int iRow,int jRow))[9]
{
    CopyMatrix(MatrixInput,MatrixInter,Row,Col);
    for(col=1;col<=Col;col++)
    {
        RowTemp[iRow][col]=MatrixInter[iRow][col];
        MatrixInter[iRow][col]=MatrixInter[jRow][col];
        MatrixInter[jRow][col]=RowTemp[iRow][col];
    }
    return MatrixInter;
}
//Pivote Downward
double MatrixDown[9][9];
double (*(PivoteDown)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixDown,Row,Col);
    Coefficient=MatrixDown[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixDown[tRow][col]/=Coefficient;
    if(tRow<Row)
        for(row=tRow+1;row<=Row;row++)
        {
            Coefficient=MatrixDown[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixDown[row][col]-=Coefficient*MatrixDown[tRow][col];
        }
return MatrixDown;
}
//Pivote Upward
double MatrixUp[9][9];
double (*(PivoteUp)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixUp,Row,Col);
    Coefficient=MatrixUp[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixUp[tRow][col]/=Coefficient;
    if(tRow>1)
        for(row=tRow-1;row>=1;row--)
        {
            Coefficient=MatrixUp[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixUp[row][col]-=Coefficient*MatrixUp[tRow][col];
        }
    return MatrixUp;
}
//Pivote in Determinant
double MatrixPiv[9][9];
double (*(Pivote)(double MatrixInput[9][9],int Dim,int pTarget))[9]
{
    CopyMatrix(MatrixInput,MatrixPiv,Dim,Dim);
    for(row=pTarget+1;row<=Dim;row++)
    {
        Coefficient=MatrixPiv[row][pTarget]/MatrixPiv[pTarget][pTarget];
        for(col=1;col<=Dim;col++)
        {
            MatrixPiv[row][col]-=Coefficient*MatrixPiv[pTarget][col];
        }
    }
    return MatrixPiv;
}
//Determinant of Square Matrix
int dCounter,dRow;
double Det;
double MatrixDet[9][9];
double Determinant(double MatrixInput[9][9],int Dim)
{
    CopyMatrix(MatrixInput,MatrixDet,Dim,Dim);
    Det=1.0;
    if(Dim>1)
    {
        for(dRow=1;dRow<Dim;dRow++)
        {
            dCounter=dRow;
            while((MatrixDet[dRow][dRow]==0.0)&(dCounter<=Dim))
            {
                dCounter++;
                Det*=-1.0;
                CopyPointer(InterchangeRow(MatrixDet,Dim,Dim,dRow,dCounter),MatrixDet,Dim,Dim);
            }
            if(MatrixDet[dRow][dRow]==0)
            {
                Det=0.0;
                break;
            }
            else
            {
                Det*=MatrixDet[dRow][dRow];
                CopyPointer(Pivote(MatrixDet,Dim,dRow),MatrixDet,Dim,Dim);
            }
        }
        Det*=MatrixDet[Dim][Dim];
    }
    else Det=MatrixDet[1][1];
    return Det;
}
//Matrix Identity
double MatrixIdent[9][9];
double (*(Identity)(int Dim))[9]
{
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            if(row==col)
                MatrixIdent[row][col]=1.0;
            else
                MatrixIdent[row][col]=0.0;
    return MatrixIdent;
}
//Join Matrix to be Augmented Matrix
double MatrixJoin[9][9];
double (*(JoinMatrix)(double MatrixA[9][9],double MatrixB[9][9],int Row,int ColA,int ColB))[9]
{
    Col=ColA+ColB;
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            if(col<=ColA)
                MatrixJoin[row][col]=MatrixA[row][col];
            else
                MatrixJoin[row][col]=MatrixB[row][col-ColA];
    return MatrixJoin;
}
//Inverse of Matrix
double (*Pointer)[9];
double IdentMatrix[9][9];
int Counter;
double MatrixAug[9][9];
double MatrixInv[9][9];
double (*(Inverse)(double MatrixInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+Dim;
    Pointer=Identity(Dim);
    CopyPointer(Pointer,IdentMatrix,Dim,Dim);
    Pointer=JoinMatrix(MatrixInput,IdentMatrix,Dim,Dim,Dim);
    CopyPointer(Pointer,MatrixAug,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            MatrixInv[row][col]=MatrixAug[row][col+Dim];
    return MatrixInv;
}
//Gauss-Jordan Elemination
double MatrixGJ[9][9];
double VectorGJ[9][9];
double (*(GaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+1;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,1);
    CopyPointer(Pointer,MatrixGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=1;col++)
            VectorGJ[row][col]=MatrixGJ[row][col+Dim];
    return VectorGJ;
}
//Generalized Gauss-Jordan Elemination
double MatrixGGJ[9][9];
double VectorGGJ[9][9];
double (*(GeneralizedGaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim,int vCol))[9]
{
    Row=Dim;
    Col=Dim+vCol;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,vCol);
    CopyPointer(Pointer,MatrixGGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(row=1;row<=Row;row++)
        for(col=1;col<=vCol;col++)
            VectorGGJ[row][col]=MatrixGGJ[row][col+Dim];
    return VectorGGJ;
}
//Matrix Sparse, Three Diagonal Non-Zero Elements
double MatrixSpa[9][9];
double (*(Sparse)(int Dimension,double FirstElement,double SecondElement,double ThirdElement))[9]
{
    MatrixSpa[1][1]=SecondElement;
    MatrixSpa[1][2]=ThirdElement;
    MatrixSpa[Dimension][Dimension-1]=FirstElement;
    MatrixSpa[Dimension][Dimension]=SecondElement;
    for(int Counter=2;Counter<Dimension;Counter++)
    {
        MatrixSpa[Counter][Counter-1]=FirstElement;
        MatrixSpa[Counter][Counter]=SecondElement;
        MatrixSpa[Counter][Counter+1]=ThirdElement;
    }
    return MatrixSpa;
}

Copy and save the above code as Matrix.h then try the following code:

#include<iostream>
#include<conio.h>
#include"Matrix.h"
int Dim;
double Matrix[9][9];
int main()
{
    cout<<"Enter your matrix dimension: ";
    cin>>Dim;
    Input(Matrix,Dim,Dim);
    cout<<"Your matrix:"<<endl;
    Output(Matrix,Dim,Dim);
    cout<<"The inverse:"<<endl;
    Output(Inverse(Matrix,Dim),Dim,Dim);
    getch();
}

Answer 5

Why don't you try to code it yourself? Take it as a challenge. :)

For a 3×3 matrix

alt text http://mathworld.wolfram.com/images/equations/MatrixInverse/NumberedEquation3.gif

the matrix inverse is

alt text http://mathworld.wolfram.com/images/equations/MatrixInverse/NumberedEquation4.gif

I'm assuming you know what the determinant of a matrix |A| is.

Images (c) Wolfram|Alpha and mathworld.wolfram (06-11-09, 22.06)

Answer 6

Don't try to do this yourself if you're serious about getting edge cases right. So while they many naive/simple methods are theoretically exact, they can have nasty numerical behavior for nearly singular matrices. In particular you can get cancelation/round-off errors that cause you to get arbitrarily bad results.

A "correct" way is Gaussian elimination with row and column pivoting so that you're always dividing by the largest remaining numerical value. (This is also stable for NxN matrices.). Note that row pivoting alone doesn't catch all the bad cases.

However IMO implementing this right and fast is not worth your time - use a well tested library and there are a heap of header only ones.