1

I am working on a c++ codebase right now which uses a matrix library to calculate various things. One of those things is calculating the inverse of a matrix. It uses gauss elimation to achieve that. But the result is very inaccurate. So much so that multiplying the inverse matrix with the original matrix isn't even close the the identity matrix.

Here is the code that is used to calculate the inverse, the matrix is templated on a numerical type and the rows and columns:

/// \brief Take the inverse of the matrix.
/// \return A new matrix which is the inverse of the current one.
matrix<T, M, M> inverse() const
{
    static_assert(M == N, "Inverse matrix is only defined for square matrices.");

    // augmented the current matrix with the identiy matrix.
    auto augmented = this->augment(matrix<T, M, M>::get_identity());

    for (std::size_t i = 0; i < M; i++)
    {
        // divide the current row by the diagonal element.
        auto divisor = augmented[i][i];
        for (std::size_t j = 0; j < 2 * M; j++)
        {
            augmented[i][j] /= divisor;
        }

        // For each element in the column of the diagonal element that is currently selected
        // set all element in that column to 0 except the diagonal element by using the currently selected row diagonal element.

        for (std::size_t j = 0; j < M; j++)
        {
            if (i == j)
            {
                continue;
            }

            auto multiplier = augmented[j][i];
            for (std::size_t k = 0; k < 2 * M; k++)
            {
                augmented[j][k] -= multiplier * augmented[i][k];
            }
        }
    }

    // Slice of the the new identity matrix on the left side.
    return augmented.template slice<0, M, M, M>();
}

Now I have made a unit test which test if the inverse is correct using pre computed values. I try two matrices one 3x3 and one 4x4. I used this website to compute the inverse: https://matrix.reshish.com/ and they do match to a certain degree. since the unit test does succeed. But once I calculate the original matrix * the inverse nothing even resembling an identity matrix is achieved. See the comment in the code below.

BOOST_AUTO_TEST_CASE(matrix_inverse)
{
    auto m1 = matrix<double, 3, 3>({
        {7, 8, 9},
        {10, 11, 12},
        {13, 14, 15}
    });

    auto inverse_result1 = matrix<double,3, 3>({
        {264917625139441.28, -529835250278885.3, 264917625139443.47},
        {-529835250278883.75, 1059670500557768, -529835250278884.1},
        {264917625139442.4, -529835250278882.94, 264917625139440.94}
    });

    auto m2 = matrix<double, 4, 4>({
        {7, 8, 9, 23},
        {10, 11, 12, 81},
        {13, 14, 15, 11},
        {1, 73, 42, 65}
    });

    auto inverse_result2 = matrix<double, 4, 4>({
        {-0.928094660194201, 0.21541262135922956, 0.4117111650485529, -0.009708737864078209},
        {-0.9641231796116679, 0.20979975728155775, 0.3562651699029188, 0.019417475728154842},
        {1.7099261731391882, -0.39396237864078376, -0.6169346682848 , -0.009708737864076772 },
        {-0.007812499999999244, 0.01562499999999983, -0.007812500000000278, 0}
    });


    // std::cout << (m1.inverse() * m1) << std::endl;
    // results in
    // 0.500000000     1.000000000     -0.500000000
    // 1.000000000     0.000000000     0.500000000
    // 0.500000000     -1.000000000    1.000000000
    // std::cout << (m2.inverse() * m2) << std::endl; 
    // results in
    // 0.396541262     -0.646237864    -0.689016990    -2.162317961
    // 1.206917476     2.292475728     1.378033981     3.324635922
    // -0.884708738    -0.958737864    -0.032766990    -3.756067961
    // -0.000000000    -0.000000000    -0.000000000    1.000000000

    BOOST_REQUIRE_MESSAGE(
        m1.inverse().fuzzy_equal(inverse_result1, 0.1) == true,
        "3x3 inverse is not the expected result."
    );

    BOOST_REQUIRE_MESSAGE(
        m2.inverse().fuzzy_equal(inverse_result2, 0.1) == true,
        "4x4 inverse is not the expected result."
    );
}

I am at my wits end. I am by no means a specialist on matrix math since I had to learn it all on the job but this really is stumping me.

The complete code matrix class is available at: https://codeshare.io/johnsmith

Line 404 is where the inverse function is located.

Any help is appreciated.

  • 2
    Some matrices don't have an inverse matrix. Like the first 3x3 example. For them, the program will produce nonsense results. – anatolyg Sep 25 '17 at 02:21
  • 2
    Make sure the determinant of the matrix is non-zero (i.e., the matrix is non-singular and can be inverted). Even if the determinant it's non-zero but very close to zero, you have to be careful (read about condition number), see e.g. https://math.stackexchange.com/questions/1622610/when-is-inverting-a-matrix-numerically-unstable – vsoftco Sep 25 '17 at 02:28
  • @vsoftco in both cases the determinant is non-zero though. –  Sep 25 '17 at 02:35
  • @anatolyg how can I know if there is an inverse? –  Sep 25 '17 at 02:35
  • @JohnSmith Is it close to zero though? In that case you should expect in-accuracies. PS: the matrix is non-invertible if the determinant is zero. It is invertible otherwise (assuming a square matrix). – vsoftco Sep 25 '17 at 02:36
  • Never heard about the side you are using to calculate the inverse, but mine says your matrix is singular http://m.wolframalpha.com/input/?i=invert+%7B7%2C+8%2C+9%7D%2C+++++++++%7B10%2C+11%2C+12%7D%2C+++++++++%7B13%2C+14%2C+15%7D – ead Sep 25 '17 at 04:09
  • Read about matrix condition numbers. https://en.wikipedia.org/wiki/Condition_number – doug Sep 25 '17 at 04:12
  • 2
    There is another issue: when using this method to invert a matrix, it's possible for the diagonal element (that will be used as a divisor) of a row to be zero, in which case that row needs to be swapped with a later row that doesn't have a zero in that column. The code will need to be updated to deal with this case. – rcgldr Sep 25 '17 at 05:24
  • Re @rcgldr comment -- this is called pivoting. In general, matrix solvers will need some form of pivoting to be numerically stable, even on well-conditioned matrices: https://en.wikipedia.org/wiki/Pivot_element – comingstorm Sep 25 '17 at 06:19
  • on top of all the comments i simply need to add: for Gauss **GEM** you need to carefully select which row/column to use for which diagonal element and in what order otherwise your result might stop on dead end for some matrices. If you need more accuracy then compute your matrix by sub-determinants. It is straight forward and uses usually bit less operations per matrix element leading to better accuracy on FPU despite the division. See [Understanding 4x4 homogenous transform matrices](https://stackoverflow.com/a/28084380/2521214) and look for `matrix_inv` implementation – Spektre Sep 25 '17 at 06:38
  • In case your matrix is orthogonal (which is most likely not the case) you can also use transpose (see the full pseudoinverse in one of the last links there for homogenuous matrices) which is fastest simplest and most precise from all the methods I know of. – Spektre Sep 25 '17 at 06:45

1 Answers1

4

As already established in the comments the matrix of interest is singular and thus there is no inverse.

Great, your testing found already the first issue in the code - this case isn't handled properly, no error is raised.

The bigger problem is, that this is not easy to detect: If there where no errors due to rounding errors, it would be a cake of piece - just test that divisor isn't 0! But there are rounding errors in floating operations, so divisor will be a very small nonzero number.

And there is no way to tell, whether this nonzero value due to rounding errors or to the fact that the matrix is near singular (but not singular). However, if matrix is near singular it has a poor condition and thus the results cannot be trusted anyway.

So ideally, the algorithm should not only calculate the inverse, but also (estimate) the condition of the original matrix, so the caller can react upon a bad condition.

Probably it is wise to use well-known and well-tested libraries for this kind of calculation - there is a lot to be considered and what can be done wrong.

ead
  • 32,758
  • 6
  • 90
  • 153