Inverse of matrix inaccurate

Question

I'm trying to write tests to validate inverse of matrix, but test fails when I use FLT_EPSILON to compare results.

Comparing function looks like this: test_assert_mat4_eq

What I'm trying to do is:

A = random matrix (4x4 float matrix)
B = inv(A)
C = inv(B)

assert(A == C) <---- fails

First let me explain how I calculate matrix inverse; mat4 is 4x4 float matrix, if SIMD is enabled then inverse of matrix will be computed via SIMD instructions (SSE2 and AVX).

You can see matrix inverse code at glm_mat4_inv if SSE2 is enabled then the inverse is calculated via glm_mat4_inv_sse2 also there is glm_mat4_inv_precise_sse2 version to get more precision by avoid _mm_rcp_ps instruction. I used the second (glm_mat4_inv_precise_sse2) version to test.

assert_true(fabsf(m1[i][j] - m2[i][j]) <= 0.0001);

passes on my macbook but it still fails on linux. 

assert_true(fabsf(m1[i][j] - m2[i][j]) <= FLT_EPSILON); this even not passes on macos. Maybe comparing with 0.001 will also work on linux but precision is too low.

I also created an issue on glm repo (https://github.com/g-truc/glm/issues/700) because this issue also valid for glm.

What is wrong with this? Why the precision is too low? Is this ok? Should I leave it like this (by removing the test or by changing precision)?

NOTE: Random matrix is generating with test_rand_mat4 function. But I'm using it just for generating a matrix. I'm not using any random matrix in anywhere, all matrices are affine transformation matrices, maybe I should use affine transforms (which is the main purpose) for comparing

Answer 1

1. A 4⨉4 matrix is small enough that you could instrument the code to print each floating-point number involved after each step of calculation and then compare the macOS and Linux results to see where they differ. (Use the %a format in C’s printf to print hexadecimal floating-point or use something like %.30g to print the entire value in decimal.)

2. Telling us the difference exceeds .0001 is meaningless since we do not know what magnitude your data is. Commonly, numerical errors will be roughly proportional to the magnitudes of some of the numbers involved.

3. If the “same” operations are performed on macOS and on Linux running on hardware that uses IEEE-754 floating-point, then identical results should be achieved. So any differences are likely because you are not using the same operations. Causes of such differences may include that the two systems are using different source code (e.g., because one uses SIMD because SSE is enabled, and the other uses scalar code) or the compilers are compiling code differently. Override the SSE decision temporarily so that you are testing the same source code on both systems. Once that is debugged, test the SSE code on both systems, if possible. After that, compare the SSE code to the non-SSE code.

4. Instead of testing with random data, you should start with simple test cases and work up to more involved test cases. Start with the identity matrix. Simple cases would help debug the fundamental logic without involving significant floating-point errors. Then modify elements to generate more complicated cases. When generating matrices for testing, avoid ill-conditioned matrices. I am not experienced in generating matrices for testing matrix inverse, so you will have to investigate how to modify a matrix to improve its condition number, or someone else might suggest something.

5. Do not test matrix inverse by calling inv twice. This misses simple errors such as a copy-and-paste error that results in a routine named inv actually being a matrix copy or a negate. Instead, use a known-good reference implementation of matrix inverse or use known test cases (the inverse of the identity matrix is the identity matrix, and others can be constructed) or use other properties of the inverse (multiplying the inverse by the original matrix should produce the identity matrix). Use double for the test code. The number of test cases needed to find errors in a 4⨉4 matrix inversion is not large enough that performance matters.

Answer 2

I double checked test codes and there was a copy-paste error in test codes. The _mm_rcp_ps instruction is used in matrix inverse function to improve performance. The problem was that I expected same precision in tests for both _mm_rcp_ps version and non-_mm_rcp_ps version. I've fixed tests.

Now assert_true(fabsf(m1[i][j] - m2[i][j]) <= 0.0000009f); passes for both on macOS and Linux for precise/accurate version

and assert_true(fabsf(m1[i][j] - m2[i][j]) <= 0.0009f) passes for both on macOS and Linux for fast (_mm_rcp_ps) version.

Same tests also passes on travis-ci.

On the other side, it wouldn't work because I used totally random matrices, even random scales wouldn't work because of being ill-conditioned matrices (thanks @Eric Postpischil and @geza)