5

The title explains it already. If I need to find an inverse of a matrix, is there any reason I should use A\eye(size(A)) instead of inv(A)? And before you ask: Yes, I really need the inverse, not only for calculations.

PS: 

isequal(inv(A), A\eye(size(A)))
ans =
 0

So which one is more precise?

UPDATE: This question was closed as it appeard to be a duplicate of the question "why is inv in MATLAB so slow and inaccurate". This question here differs significantly by not addressing the speed, nor the accuarcy of the function inv but the difference of inv and .\eye to calculate the true inverse of a matrix.

Community
  • 1
  • 1
rst
  • 2,510
  • 4
  • 21
  • 47
  • The inv seems to be a bit faster for me on a rand(4000) matrix but I have no explanation why. The difference in infinite norm is 1e-10. It is really an interesting question. – Zoltán Csáti Jan 27 '16 at 10:15
  • Just out of curiosity, what do you need the inverse for? – Rody Oldenhuis Jan 27 '16 at 11:35
  • 1
    More specifically, why do you need the inverse to be precise? cc: @RodyOldenhuis – Stewie Griffin Jan 27 '16 at 11:36
  • 1
    The inverse in a computer is generally a bad idea! ` \ `<- this does very advanced math stuff in order to have a better solution – Ander Biguri Jan 27 '16 at 12:11
  • @RodyOldenhuis I won't to into details but we have a optimization scheme, where to calculate and benchmark a sensitivity, we need to have the inverse of a complex matrix. I already managed to overcome the need of calculating the actual inverse, however, the question I thought was interesting. – rst Jan 27 '16 at 12:17
  • There are some good answers, which one to accept? – rst Jan 27 '16 at 12:18

3 Answers3

8

Let's disregard performance (speed) and best practice for a bit.

eps(n) is a command that returns the distance to the next larger double precision number from n in MATLAB. So, eps(1) = 2.2204e-16 means that the first number after 1 is 1 + 2.2204e-16. Similarly, eps(3000) = 4.5475e-13. Now, let's look at the precision of you calculations:

n = 100;
A = rand(n);
inv_A_1 = inv(A);
inv_A_2 = A \ eye(n);

max(max(abs(inv_A_1-inv_A_2)))
ans =
   1.6431e-14

eps(127) = 1.4211e-14
eps(128) = 2.8422e-14

For integers, the largest number you can use that has an accuracy higher than the max difference between your two matrices is 127.

Now, let's check how the accuracy when we try to recreate the identity matrix from the two inverse matrices.

error_1 = max(max(abs((A\eye(size(A))*A) - eye(size(A)))))
error_1 =
   3.1114e-14
error_2 = max(max(abs((inv(A)*A) - eye(size(A)))))
error_2 =
   2.3176e-14

The highest integer with a higher accuracy than the maximum difference between the two approaches is 255.

In summary, inv(A) is more accurate, but once you start using the inverse matrices, they are for all intended purposes identical.

Now, let's have a look at the performance of the two approaches:

n = fix(logspace(1,3,40));
for i = 1:numel(n)
A = rand(round(n(i)));
t1(i) = timeit(@()inv(A));
t2(i) = timeit(@()A\eye(n(i)));
end
loglog(n,[t1;t2])

enter image description here

It appears that which of the two approaches is fastest is dependent on the matrix size. For instance, using inv is slower for n = 255, but faster for n = 256.

In summary, choose approach based on what's important to you. For most intended purposes, the two approaches are identical.

Note that svd and pinv may be of interest if you're working with badly scaled matrices. If it's really really important you should consider the Symbolic toolbox.

I know you said that you "actually need the inverse", but I can't let this go unsaid: Using inv(A)*b is never the best approach for solving a linear equation! I won't explain further as I think you know this already.

Stewie Griffin
  • 14,889
  • 11
  • 39
  • 70
4

If you need the inverse, you should use inv.

The inverse is calculated via LU decomposition, whereas the backslash operator mldivide calculates the solution to your linear system using different methods depending on the properties of your matrix A (see https://scicomp.stackexchange.com/a/1004), which can yield less accurate results for the inverse.

It should be noted that if you want to solve a linear system, the calculation is likely going to be much faster and more accurate using mldivide(\). The MATLAB documentation of inv is basically one big warning not to use inv to solve linear systems.

Community
  • 1
  • 1
dasdingonesin
  • 1,347
  • 1
  • 10
  • 16
2

Just a way trying to check this, not sure if it's completely helpful though: multiply your inverse matrix result back with it's original version and check the deviation from the identity matrix:

A = rand( 111 );
E1 = abs( (A\eye(size(A) ) * A ) - eye( size(A) ) );
E2 = abs( ( inv(A) * A )         - eye( size(A) ) );
mean(E1(:))
mean(E2(:))

inv seems to be more accurate as I would have expected. Maybe somebody can re-evaluate this. ;)

Matthias W.
  • 1,039
  • 1
  • 11
  • 23
  • Can't it be a cancellation error, because both inv(A)*A and A\eye(size(A))*A is close to eye(size(A))? However, the cancellation error would be present for both cases, so inv seems to be more accurate. A floating point arithmetic exercise would help, but I am not proficient in that. – Zoltán Csáti Jan 27 '16 at 10:22
  • I get the same result, E2 having slightly less mean value. – rst Jan 27 '16 at 12:04