Accuracy of matrix inverse

[martinv]

Hello,

I use dgetrf and dgetri (called from C++) to invert the following matrix:

Code: Select all

A = [
1 -1 -1 1 1 1
1 1 -1 1 -1 1
1 0 2 0 0 3
1 0 -1 -0.5 -0 1
1 0.5 0.5 0.25 0.75 -0.5
1 -0.5 0.5 0.25 -0.75 -0.5]


The inverse that I obtained is:

Code: Select all

inv(A) = [
-2.77556e-17 5.55112e-17 -2.77556e-17 0.333333 0.333333 0.333333
-0.3 0.3 0 0 0.4 -0.4
-0.1 -0.1 0.2 -0.266667 0.133333 0.133333
0.333333 0.333333 0 -0.666667 0 0
0.2 -0.2 0 0 0.4 -0.4
0.0666667 0.0666667 0.2 0.0666667 -0.2 -0.2 ]


What disturbs me are the near-zero values (effect of finite precision artihmetic, I know). The problem is that I don't see the same when running Matlab:

Code: Select all

inv(A)
ans =
         0   -0.0000         0    0.3333    0.3333    0.3333
   -0.3000    0.3000         0         0    0.4000   -0.4000
   -0.1000   -0.1000    0.2000   -0.2667    0.1333    0.1333
    0.3333    0.3333         0   -0.6667         0         0
    0.2000   -0.2000         0         0    0.4000   -0.4000
    0.0667    0.0667    0.2000    0.0667   -0.2000   -0.2000


Is there anything I can do to increase the accuracy of the matrix inverse? I am using the reference Lapack implementation installed through my Linux distribution.

Thank you very much.

Best regards,
Martin Vymazal

[Julien Langou]

Hi Martin,

You are fine, no idea why Matlab gives a ZERO and LAPACK gives you a 2e-17. This is not a problem of increasing the precision of LAPACK. MATLAB is using LAPACK in any case. Why does MATLAB differs slightly from LAPACK on this instance, no idea. Maybe you and Matlab are using different BLAS, etc. Your C++ call to LAPACK is fine and LAPACK is accurate enough. In other world, if MATLAB has one ZERO correct and LAPACK does not, this is essentially luck. It could as well be the reverse next time.

J.

[martinv]

Hello Julien,

I'm sorry, but the fault is actually on my side. I checked again and it seems that some of the zero entries in A were not exactly equal to zero. I didn't print them on the screen with sufficient precision, so they only appeared to be zero. In the end, Lapack did a correct job ...

Sorry for the misleading post ...

Best regards,

Martin Vymazal

[Julien Langou]

No worries. I would check the entry (1,1) of inv(A) with: B=inv(A); B(1,1)==0 and this would say whether the entry (1,1) of inv(A) is 0 or not. No big deal. Thanks for posting. J.