Trust me, you do not want to see the formulas for higher-order determinants. But as you might surmise, there exist recursive algorithms to compute them from lower-order determinants. Writing the computer code is easier than writing the math equation.

The important thing about determinants and singular matrices is that, just as the CPU catches an attempt to divide by zero, we must catch an attempt to invert a singular matrix. In practice, this means we must calculate the determinant as we calculate the inverse. If it ever becomes zero, we must call a halt and report failure.

When a CPU detects an attempt to divide by zero, it issues a hardware exception. In a perfect world, we should be able to throw a software exception for the singular matrix, as well. However, the C++ exception mechanism is complicated and usually not allowed in embedded software. We also can't just send an error message to the console. In general, embedded systems don't have consoles. So instead, I'm returning a flag of sorts: the determinant value. If it's zero, the inversion failed.

mInv, the matrix inversion algorithm I'm providing in SimpleMat.cpp, is an oldie but goodie. I cribbed it from the (uncopyrighted) IBM Scientific Subroutine Package, ca. 1960. I've discussed this algorithm at length in earlier issues of Embedded Systems Design, so I won't go into it in detail--perhaps later. The algorithm uses Gauss-Jordan elimination with full pivoting. This makes it something of an oddity in the world of numerical analysis.

About pivoting: as the Gauss-Jordan process unfolds, the key elements turn out to be those down the diagonal of the matrix. We end up dividing by those elements, so needless to say, zeros on the diagonal are bad. In fact, near-zero values are also bad. To get the best numerical accuracy, we'd prefer to have large (positive or negative) elements there. The term "pivoting" refers to a process of seeking out the larger elements and moving them to the diagonal, by row, and perhaps column, transformations.

I have seen, in embedded systems, homebrewed matrix inversion software that did no pivoting at all. This is both a terrible idea and a perfect example of what can happen when embedded programmers go bad. It's perfectly possible for a matrix to have a zero on the diagonal but still be non-singular. By performing no pivoting, the programmers risk having an algorithm fail when it doesn't need to.

Most implementations of the Gauss-Jordan algorithm use partial pivoting: they seek out the largest element but only in a given column. Full pivoting looks in the remaining elements of both rows and columns, thereby assuring the largest possible element at each step. Full pivoting takes a little longer. When you're moving elements to the diagonal, you're shuffling the matrix, much like a Rubik's cube. As you do so, you must record each shuffle, so that at the end of the process, you can put it back into order. Because full pivoting uses both row and column transformations, the shuffle/ unshuffle process takes longer and uses a bit more memory.

The upside is, you can get larger elements down that critical diagonal, and therefore the highest accuracy possible. The IBM programmers thought this was worth the extra work, and I agree. Which is why I've been using the algorithm for lo these 40+ years.

I should mention that Gauss-Jordan is not the only way to invert a matrix. More modern methods involve factoring the matrix into two parts. The LU decomposition, for example, involves factoring the matrix into two triangular matrices. This process gives a little better performance if you're really seeking the matrix inverse, as you do to solve Equation 7 via Equation 8. However, if you have to solve many such linear equations, involving the same matrix A, LU decomposition is much faster.

Several other decomposition algorithms exist. Each has its advantages. However, if you want the matrix inverse explicitly, and your matrix is relatively small (say less than 20x20), you'll find mInv hard to beat.