0

I am calculating the inverse of a square matrix using different libraries via Cholesky factorization. However my results are not as I was expecting. I am not an expert in maths, but I was expecting to get a closer result.

I am using MLK, magma and CULA libraries to calculate the inverse of a matrix in CPU and GPUs. After doing the calculation those libraries I've noticed the results always differ in one element. Say I want to calculate the inverse of A= [0.237306,0.000458;0.000458,0.238497]:

A[0] = 0.237306 
A[1] = 0.000458 
A[2] = 0.000458 
A[3] = 0.238497

The result I obtain is:

inv(A)[0] = 4.213983 
inv(A)[1] = -0.008092 
inv(A)[2] = 0.000458 
inv(A)[3] = 4.192946

However, the correct result should be 

   4.2139841  -0.0080924
  -0.0080924   4.1929404

As you can see, inv(A)[3] is different, although the rest of them are fine. Is that how Cholesky Inversion should work? Is this a correct/approximate result or am I doing something wrong here?

Amro
  • 123,847
  • 25
  • 243
  • 454
Manolete
  • 3,431
  • 7
  • 54
  • 92
  • Either provide information on used libraries, or try http://math.stackexchange.com/ – ipavlic Jul 20 '12 at 12:02
  • 7
    Are you sure that the values are all supposed to be usable ? Typically symmetric matrix libraries use only as output one side of the matrix, relying on you to complete by symmetry the other side. Please post the outline of the code leading to the result. – Alexandre C. Jul 20 '12 at 12:07

2 Answers2

3

I understand what is going on. This libraries calculate the inverse of the matrix, upper or lower as Alexandre C pointed out. So depending on the argument you pass to the library (Upper or Lower), it calculates the inverse of the upper or lower side of the matrix. I thought you could calculate the complete inverse matrix, but apparently it is not possible this way.

Manolete
  • 3,431
  • 7
  • 54
  • 92
1

It's probably due to the IEEE spec for floating point numbers and the imprecision that's built into working with them.

You don't say if you've specified 32-bit or 64-bit floating point numbers. Try increasing your precision.

Another point to make is that people don't usually calculate the inverse of a matrix because it's prone to suffer from round off errors. A better choice is LU decomposition and forward back substitution.

Where did you get the "correct" answer from? How are you judging the correctness of your result?

I know it's not the best tool for numerical methods, but if I take your matrix:

0.237306 0.000458 0.000458 0.238497

and run it through MINVERSE, I get this answer:

4.213984 -0.00809 -0.00899 4.19294

I can check its validity by multiplying the original and the inverse to see if I get an identify matrix back. Here's what I got out of MMULT:

1 2.1684e-19 -2.168e-19 1

That's what a numerical methods person would call an identity matrix.

So I say the answer I get is correct. You should check yours in the same way, and read What Every Computer Scientist Should Know About Floating Point Arithmetic.

duffymo
  • 305,152
  • 44
  • 369
  • 561