Calculating the Inverse of a 3 by 3 Matrix

amalfiCoast

Hi Everyone,

I'm having a problem calculating the inverse of a matrix. My program is not yet complete but my approach is to first calculate the determinant of a matrix, then the adjoint of the matrix, then multiply the determinant and the adjoint to get the inverse. The error I'm getting is within the matrixInverse function at

adj_p[3][3] = {

"error: expected expression before '{' token". Am I not allowed to multiply terms inside an array?

Also, this code as well as my other matrixMultiplication program (posted a month ago) seem to be hard-coded for 3 by 3 matrices only. What if I have a situation where I need to manipulate a 2 by 2 matrix, for example? Do I need to make other functions for other dimensioned matrices or is there a standard approach that works for all matrix sizes? In Matlab, for example, you don't need to specify or hard-code the matrix dimension. By the way, I'm developing these functions because I want to implement an Extended Kalman Filter on hardware. We did a project in school designing and simulating a UAV autopilot in Matlab and I've always wanted to be able to implement EKFs in hardware.

/* Take the inverse of a 3 by 3 matrix */
#include "simpletools.h"

float p[3][3] = {       // Define matrix p
  {1, 1, 2},
  {3, 4, 5},
  {6, 7, 8}
};

float r[3][3] = {       // Zero matrix r
  {0, 0, 0},
  {0, 0, 0},
  {0, 0, 0}
};

float det_p[1];         // Define determinant of p (scalar)

float adj_p[3][3] = {       // Zero matrix adj_p
  {0, 0, 0},
  {0, 0, 0},
  {0, 0, 0}
};

float matrixInverse(float det_p[1], float adj_p[3][3], float p[3][3], float r[3][3]);        // Declare matrixInverse function

int main(){
  
  matrixInverse(det_p, adj_p, p, r);        // Call matrixInverse function

  // Print matrix adj_p 
  print("adj_p =\n"); 
  for(int i = 0; i < 3; i++){
    for(int k = 0; k < 3; k++){
      print("%f ", r[3][3]);
    }   
    print("\n");
  }
  // Print determinant of p
  print("\ndet_p = %f\n", det_p[1]);
}  

// matrixInverse function - compute the inverse of a 3 by 3 matrix
float matrixInverse(float det_p[1], float adj_p[3][3], float p[3][3], float r[3][3]){
  
  det_p[1] = p[0][0]*p[1][1]*p[2][2]+p[0][1]*p[1][2]*p[2][0]+p[0][2]*p[1][0]*p[2][1]    // Calculate determinant of p
            -p[0][2]*p[1][1]*p[2][0]-p[0][1]*p[1][0]*p[2][2]-p[0][0]*p[1][2]*p[2][1];
  
  adj_p[3][3] = {
    {p[1][1]*p[2][2]-p[1][2]*p[2][1], p[0][2]*p[2][1]-p[0][1]*p[2][2], p[0][1]*p[1][2]-p[1][1]*p[0][2]}
    {p[2][0]*p[1][2]-p[1][0]*p[2][2], p[0][0]*p[2][2]-p[2][0]*p[0][2], p[1][0]*p[0][2]-p[0][0]*p[1][2]},
    {p[1][0]*p[2][1]-p[2][0]*p[1][1], p[2][0]*p[0][1]-p[0][0]*p[2][1], p[0][0]*p[1][1]-p[0][1]*p[1][0]}
  };    
}

Thank you,
David

Roy Eltham

  adj_p[3][3] = {
    {p[1][1]*p[2][2]-p[1][2]*p[2][1], p[0][2]*p[2][1]-p[0][1]*p[2][2], p[0][1]*p[1][2]-p[1][1]*p[0][2]}
    {p[2][0]*p[1][2]-p[1][0]*p[2][2], p[0][0]*p[2][2]-p[2][0]*p[0][2], p[1][0]*p[0][2]-p[0][0]*p[1][2]},
    {p[1][0]*p[2][1]-p[2][0]*p[1][1], p[2][0]*p[0][1]-p[0][0]*p[2][1], p[0][0]*p[1][1]-p[0][1]*p[1][0]}
  };

This syntax is only valid for initializing, not for general assignment. You need to do this as the individual assignments of the adj_p elements.

 adj_p[0][0] = p[1][1]*p[2][2]-p[1][2]*p[2][1];
 adj_p[0][1] = p[0][2]*p[2][1]-p[0][1]*p[2][2];
 ...
 adj_p[2][2] = p[0][0]*p[1][1]-p[0][1]*p[1][0];

Also, your det_p[1] assignment is invalid, it's accessing outside of the array. It should be det_p[0] = ... In C, array indices are 0 based, so a declaration of float x[3] is accessed with indices from 0 to 2.

amalfiCoast

Wow, thanks Roy! I was actually confused about the det_p[1] assignment so I'm glad you brought that up. It does work because it does calculate the determinant correctly so why wouldn't it be giving me an error? Is it correct in the initialization and each subsequent declaration because it describes a vector of 1 element (scalar), similar to the p[3][3] describing a 3 by 3 matrix? Should it only be det_p[0] when I print, since I want to print the term that's in the "0 position"? Thank you again so much! I sat here scratching my head for quite a while and finally decided to ask for help

David

amalfiCoast

Hi,

The code now works as it should. Please see below. I am however still confused about the det_p[1] issue that was discussed. Also, is there a way for this code to be adapted to any size matrix without having to recode the function?

Thank you again, Roy, for your help.

/* Take the inverse of a 3 by 3 matrix */
#include "simpletools.h"

float p[3][3] = {       // Define matrix p
  {1, 1, 2},
  {3, 4, 5},
  {6, 7, 8}
};

float r[3][3] = {       // Zero matrix r
  {0, 0, 0},
  {0, 0, 0},
  {0, 0, 0}
};

float det_p[1];         // Define determinant of p (scalar)

float adj_p[3][3] = {       // Zero matrix adj_p
  {0, 0, 0},
  {0, 0, 0},
  {0, 0, 0}
};

float matrixInverse(float det_p[1], float adj_p[3][3], float p[3][3], float r[3][3]);        // Declare matrixInverse function

int main(){
  
  matrixInverse(det_p, adj_p, p, r);        // Call matrixInverse function
/*
  // Print matrix adj_p 
  print("adj_p =\n"); 
  for(int i = 0; i < 3; i++){
    for(int k = 0; k < 3; k++){
      print("%f ", adj_p[i][k]);
    }   
    print("\n");
  }
  // Print determinant of p
  print("\ndet_p = %f\n", det_p[1]);
  */
  // Print matrix inverse, r 
  print("\nr =\n"); 
  for(int i = 0; i < 3; i++){
    for(int k = 0; k < 3; k++){
      print("%f ", r[i][k]);
    }   
    print("\n");
  }
}  

// matrixInverse function - compute the inverse of a 3 by 3 matrix
float matrixInverse(float det_p[1], float adj_p[3][3], float p[3][3], float r[3][3]){
  
  det_p[1] = p[0][0]*p[1][1]*p[2][2]+p[0][1]*p[1][2]*p[2][0]+p[0][2]*p[1][0]*p[2][1]    // Calculate determinant of p
            -p[0][2]*p[1][1]*p[2][0]-p[0][1]*p[1][0]*p[2][2]-p[0][0]*p[1][2]*p[2][1];
  
  adj_p[0][0] = p[1][1]*p[2][2]-p[1][2]*p[2][1];      // Calculate adjoint of p
  adj_p[0][1] = p[0][2]*p[2][1]-p[0][1]*p[2][2];
  adj_p[0][2] = p[0][1]*p[1][2]-p[1][1]*p[0][2];
  adj_p[1][0] = p[2][0]*p[1][2]-p[1][0]*p[2][2];
  adj_p[1][1] = p[0][0]*p[2][2]-p[2][0]*p[0][2];
  adj_p[1][2] = p[1][0]*p[0][2]-p[0][0]*p[1][2];
  adj_p[2][0] = p[1][0]*p[2][1]-p[2][0]*p[1][1];
  adj_p[2][1] = p[2][0]*p[0][1]-p[0][0]*p[2][1];
  adj_p[2][2] = p[0][0]*p[1][1]-p[0][1]*p[1][0];
  
  for(int i = 0; i < 3; i++){       // Calculate matrix inverse of p
    for(int k = 0; k < 3; k++){
      r[i][k] = adj_p[i][k]*1/det_p[1];
    }
  }               
}

Respectfully,
David

JasonDorie

det_p is declared as an array with 1 element, and C uses zero-based indices, so the only valid index is zero.

The compiler isn't very smart, however, and likely treating det_p as an array of "some number" of elements. When the code is compiled, det_p[1] is treated as the second index into a table of N elements, so you're accessing the 4 bytes past det_p[0]. Most likely it'll be adj_p[0][0], because that's the next thing declared, but it depends on how the compiler arranges your variables in memory. If it works, you got lucky.

In C, it's actually legal to declare an array of zero elements, like this:

  float dumbArray[0];

  dumbArray[5] = 2.0f;

I don't recommend it. :-)

AJL

amalfiCoast wrote: »

Hi,

The code now works as it should. Please see below. I am however still confused about the det_p[1] issue that was discussed.

As your determinant is a scalar, do you even need the index? I recognise that in matrix maths notation it is represented as a 1 x 1 matrix, but that is a convention that confuses the issue here. Perhaps just declare and use det_p

amalfiCoast wrote: »

Also, is there a way for this code to be adapted to any size matrix without having to recode the function?

To recode this only once for general matrix sizes, you'd need to code this matrixInverse() function with input parameters being pointers to the arrays rather than copies of the arrays, you'd need to pass an array size parameter, and your code would need to iterate over the arrays in similar way to your matrix printing code.

By identifying the patterns in the array indices used in the operations you should be able to determine algorithms to perform general size matrix operations.

An alternative is to look for code already written that serves the purpose.

Regards,
Anthony.

amalfiCoast

Thanks everyone - I figured out where my indexing error was.

Anthony: I'll work on modifying the code so that it will work for a general matrix. I'll report back when I have an update.

Thanks,
David