arm_mat_inverse_f64.c (V.1.4.5 and V.1.5.1) produces faulty results

[I3I-I]

Greetings!

I'm currently working on 66x66 matrices on a STM32F407 microcontroller during my master thesis at Technical University Munich and ran into a problem with CMSISs matrix inversion. I tried to get the inverse of a kind of tridiagonal block matrix (J_stm) but the result I got had errors. I compared the resulting inverse (J_stm_arm_inv) with the result from matlab (J_stm_matlab_inv) and a Gauß-Jordan Inversion I wrote myself.

Multiplicating the original matrix (J_stm) with its CMSIS inverse (J_stm_arm_inv) does not produce an eye matrix (J_stm_eye_arm): there are some 0.5 entries left and it is no diagonal matrix. See the picture for clarification. Using the results from matlab and my own routine produces an eye matrix with some negligible errors.

Reproducing the error with smaller matrices failed (correct results there), so I had to provide you the big matrices as csv files.

The arm_mat_inverse_f64.c version in use for generating the appended data was V.1.4.5, but V.1.5.1 did provide uncorrect results as well.

I hope someone can take the time to take a look at my files and check if she/he got's the same error.

Kind regards,

Ben H testdata_matrix_inverse.zip

[JonatanAntoni]

Hi Ben,

thanks for your detailed report. We'll take a look into it.

Best, Jonatan

[llefaucheur]

Hi Ben, thank you for the data that greatly helps debugging this issue. I tried to reproduce your issue using Matlab code below: [J_stm, n, d] = xlsread(' J_stm.csv'); % read original data [J_stm_matlab_inv, n, d] = xlsread(' J_stm_matlab_inv.csv'); % read Matlab inverse max(max(round(J_stm_arm_inv - inv(J_stm)))) % compute the inverse and find the max error

The error is about 1/100k the maximum observed in the inverted matrix (which is poor compared to the 10^15 dynamic range of double precision float) and could the sign you are pushing the algorithm to the limit of double precision accuracy. Would you be kind to reproduce this on your side ? We continue here finding the root cause of you observed. Best regards, Laurent.

[llefaucheur]

For example, the column #4 of J_stm has data with 130dB dynamic range, which pushes the arithmetics operation to the limits. STM32F407 is made of a CortexM4 with simple-precision floating-point operations. In order to do computations in double-precision, the compiler is using of a software library. Does the computation libraries are identical, if no we will need to have a look at the details.

[I3I-I]

Hi Laurent, thank you for your quick response! I tried to run the matlab skript you posted, but there were some minor errors (xslread, can't handle csv in my Matlab version (R2017a) and J_stm_arm_inv was not imported into the workspace). Here is the skript I used:

J_stm = csvread('J_stm.csv'); % read original data J_stm_matlab_inv = csvread('J_stm_matlab_inv.csv'); % read Matlab inverse J_stm_arm_inv = csvread('J_stm_arm_inv.csv'); % read stm inverse maxerr = max(max(round(J_stm_arm_inv - inv(J_stm)))); % compute the inverse and find the max error maximum = max(max(abs(J_stm_arm_inv))); % biggest value in arm inverse ratio = maxerr/maximum;

The output I get is: maxerr = 327911 maximum = 362790 ratio = 0.9039

So the error is the same as yours, but i can't reproduce the ratio of 1/100k.

The only external library I'm using for computation is math.h. The compiler is set to use the FPU. (screenshot appended).

Best regards, Ben

[I3I-I]

I might have found something here: Take a look at the original matrix J_stm and matrix J_stm_eye_arm = J_stm_arm_inv * J_stm, which should be an eye matrix. After column 29, J_stm_eye arm and J_stm are exactly the same. Column 29 is the first column in which the only entry is above the diagonal. I think your algorithmn only considers elements below the diagonal, so after column 29 there might be some unexpected behavior.

Would you be so kind to check if my thoughts might be correct?

Best regards, Ben

[I3I-I]

Hi Laurent,

did you have a chance to take a look at the matrix inversion again?

Best regards, Ben

[llefaucheur]

Hi Ben,

Sorry for the delays. Yes the problem seems in the library. I reproduced it after porting the code on a PC. The problem is also there with smaller size matrix (35x35 for example).

Best regards, Laurent

[I3I-I]

No problem, thank you for the update Laurent.

I'm happy to hear that you are working on a solution. Please tell me if I can be of any help!

Best regards, Ben

[llefaucheur]

Your problem gave us the opportunity to look at the sensitivity of the algorithm with ill-conditioned large matrix. For example please have look at this ACM paper "On the Stability of Gauss-Jordan Elimination with Pivoting" about rounding errors. Matlab inv() algorithm is using the LU decomposition which is less sensitive to errors, while our Gauss-Jordan algorithm is optimised for speed. We made other experiments using the "long double" format proposed for x86 processor, but here again the accuracy is not enough for your dataset example. I would like to propose you to close this topic and let us define an other API optimised for accuracy in F64 formats (LU decomposition or div&conquer). And we will add some comments in the documentation to link the data dynamic range I mentioned before, and the floating-point IEEE format able to cope with it.

[I3I-I]

Thank you very much for your efforts! I'll take a look at the paper you mentioned. For now I'm sticking to the algorithm I've written based on the Gauss-Jordan but I'm looking forward to the solution you propose in the future.

Kind regards, Ben