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.

Solution 1

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

For a 3×3 matrix

_{(source: wolfram.com)}

the matrix inverse is

_{(source: wolfram.com)}

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)

Solution 2

Here's a version of batty's answer, but this computes the correct inverse. batty's version computes the transpose of the inverse.

// computes the inverse of a matrix m
double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) -
             m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
             m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));

double invdet = 1 / det;

Matrix33d minv; // inverse of matrix m
minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet;
minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet;
minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet;
minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet;
minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet;
minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet;
minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet;
minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet;
minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet;

Solution 3

This piece of code computes the transposed inverse of the matrix A:

double determinant =    +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
                        -A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0))
                        +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
double invdet = 1/determinant;
result(0,0) =  (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
result(2,0) =  (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
result(1,1) =  (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
result(0,2) =  (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
result(2,2) =  (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;

Though the question stipulated non-singular matrices, you might still want to check if determinant equals zero (or very near zero) and flag it in some way to be safe.

Solution 4

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!

Solution 5

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.

View more solutions

Comments

[Not Sure]

You give no details and your question is very generic, so my answer is to use BLAS.

[davidtbernal]

When it comes to matrix inversion, 3x3, and cramer's rule are pretty detailed.

[batty]

In fairness, I added the additional detail after he complained. ;-)

[Markus Schnell]

"I'd prefer simplicity over speed" The problem you will be facing is numerical errors. Are you sure you wouldn't want to include a dependable library?

[nlucaroni]

dgetri will do the trick

[batty]

Yeah, this isn't code that I'm going to keep around. Later on I'll replace it with BLAS, but in the short term I have other priorities.

[duffymo]

Nice typesetting as well. LaTeX?

[SuPra]

As big a fan as I am of LaTeX, no. They are JPEGs generated by Wolfram|Alpha. :)

[Nick Dandoulakis]

Don't forget to check if determinant is zero :)

[sidgeon smythe]

Now try doing that without so much... mess. :) (Ideally, without using any number other than "0" and "3" in the code.)

[sidgeon smythe]

Uh, is this actually right? All the off-diagonal entries look wrong to me, but maybe I'm not aware of some interesting matrix identity. E.g., shouldn't the cofactor of a12 be -|a21 a23; a31 a33|?

[SuPra]

That is perfectly right. The way you're putting it will be right if the rows and columns are interchanged and the determinant sign is reversed.

[sourcenouveau]

What Wolfram Alpha search did you use to get them?

[Mr.Ree]

It looks like the jpegs off the wolframe math site. Try google. Yes the matrix is correct. Inverse is 1/det * minor of the TRANSPOSED matrix.

[batty]

"With all due respect" my code is crap? Nice. ;-) Maybe you could do the forward error analysis to determine some bounds on the round-off error instead of just picking an arbitrary constant cutoff for the determinant? 1e-2 is probably rather too conservative of a cutoff.

[Mr.Ree]

If the determinant is very large, 1/det (which we multiply by) is close to zero. If the determinant is very small, 1/det has divide by zero issues and becomes extremely large. Where to set the thresholds is somewhat application dependent, and perhaps best decided stochastically (probabilistically) depending on your data. Alternatively, you may want to consider libgmp/libgmpxx for improved precision.

[shoosh]

This code actually gives you the TRANSPOSE of the inverse matrix. Check out the formula at the other answer

[Oliver Turner]

"It is permitted to use and post individual, incidental results or small groups of results from Wolfram|Alpha on non-commercial websites and blogs, provided those results are properly credited to Wolfram|Alpha (see Attribution and Licensing below)." -- Where's the citation?

[AbePralle]

Concise, elegant, and fast (even if it is the transpose) - nicely done! A more general solution that uses loops and multiple function calls is both overkill and needless overhead for such a fundamental operation.

[Robinson]

The interesting thing here is that Mr Ree's code is a few more orders of magnitude more difficult to read and understand than the one given above.

[rbaleksandar]

I also have some advice that spares some unnecessary calculation. If you are also working with rotation matrices, their inverse is equal to their transpose. Rot^(-1) = Rot^t In such case the formula mentioned by @Suvesh Pratapa can be skipped, which adds a speed and memory boost if you deal with a huge amount of rotations (example: robotics).

[Matthew]

This was my try... i guess it doesn't work right... sorry. lol. I didn't try it before hand. :/

[Stefan]

Wonderful solution. I strongly prefer it to any solution that uses many function calls.

[Michael Anderson]

I believe that this method is numerically unstable for some nearly singular matrices.

[Michael Anderson]

This seems a very good choice to me - you can use Matrix33::inverse if you don't care about stability, and Matrix33::gjInverse if you do.

[Admin]

Nice if you have the time to take a challenge but frankly it doesn't answer the question: "I'm just looking for a short code snippet that'll do the trick for non-singular matrices"

[Luc Bloom]

This is the best answer imo, but it calculates some things 2 times, like m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2). Please store the 2nd part of the 3 determinant calculations in 3 temporary doubles. They're going to be used in the last part of the calculation again. If you flip m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0), you can even re-use it without the minus sign.

[Cornstalks]

I wrote this primarily for clarity, not necessarily speed. But you're absolutely right that some operations are repeated and can be avoided by temporarily caching the result and reusing it.

[moldovean]

what method did you use to compute the inverse, Say OL?

[moldovean]

you don't really need dim if you use the vector library because you can use the property .size()

[Hugo Maxwell]

The compiler should be able to optimize that. No worries.

[lightxbulb]

A quite elegant way to do it, if the matrix's rows are 3d vectors v[0],v[1],v[2], would be invdet = 1/Sum(v[0]*cross(v[1],v[2])); and temp = transpose(mat3(cross(v[1],v[2]),cross(v[2],v[0]),cross(v[0],‌​v[1]))) then the inverse matrix is det*temp. Where * is componentwise multiplication, Sum(v) := v[0]+v[1]+v[2];, cross is the cross product, mat3 is a constructor of a 3x3 matrix class that takes 3 3d vectors as arguments.

[Alan Wolfe]

I'm going to be that guy. Downvoted because this doesn't actually answer the question - the C++ code is not here.

[Tim Čas]

The link is down.

[Kenn Sebesta]

the QMatrix.h link is sadly dead, invalidating this answer. Perhaps the class could be included here as an appendix?

[Kenn Sebesta]

It's bad mojo to have a code answer which is known to be defective. @Matthew, if you can't confirm that this works, can you remove the answer?

[moldovean]

thanks for the update Kenn. It was working in 2015 :)

[Maurizio Tomasi]

This is the simplest implementation of the algorithm. However, beware that the determinants of the minors might be affected by round-off errors. This site provides a good explanation of the problem, and a nice solution: pharr.org/matt/blog/2019/11/03/difference-of-floats.html

[Константин Ван]

Hey, I can’t read those terms on the black background.

[Ivan Kovtun]

The computations are correct, but result become unstable when determinant is close to 0. I found this by comparing it with with the numpy.linalg.inv method which provides much better result. Just multiply the inverted and the original matrices and compare it with the identity matrix.