inv(matrix) takes shorter time than \ operator

A=magic(5)

A =

      17    24     1     8    15
      23     5     7    14    16
       4     6    13    20    22
      10    12    19    21     3
      11    18    25     2     9

>> tic
inv(A)
toc

ans =

     -0.0049    0.0512   -0.0354    0.0012    0.0034
      0.0431   -0.0373   -0.0046    0.0127    0.0015
     -0.0303    0.0031    0.0031    0.0031    0.0364
      0.0047   -0.0065    0.0108    0.0435   -0.0370
      0.0028    0.0050    0.0415   -0.0450    0.0111

Elapsed time is 0.006097 seconds.
>> tic
B=eye(5);
A\B
toc

ans =

     -0.0049    0.0512   -0.0354    0.0012    0.0034
      0.0431   -0.0373   -0.0046    0.0127    0.0015
     -0.0303    0.0031    0.0031    0.0031    0.0364
      0.0047   -0.0065    0.0108    0.0435   -0.0370
      0.0028    0.0050    0.0415   -0.0450    0.0111

Elapsed time is 0.008993 seconds.

Just a general question, shouldn't have been the other way round? inv(matrix) is supposed to be more computationally intensive than the \ operation.

[Jan]

This is an basic and important problem in the field of numerical linear algebra. +1

[Matt J]

Armin Schmidt's comments moved here:

Hey, I know this discussion is already a little bit older, but I have a similar question to inv and \, but in relation with the mahalanobis-distance. I can get from an other class the covariance-matrices of each cluster. Now I want to calculate the mahalanobis-distance: (x-mu)*cov^(-1)*(x-mu)', where mu is the centroid of each cluster and x the current feature-vector to measure the distance to my clusters (classes). Now to write A*inv(cov)*A' (for A=(x-mu)) would be numerically instable. what would be the best way to calculate the mahalanobis-distance in this case? Could I write something like A * cov\A' and would that be numerically more stable?

[Matt J]

@Armin,

Yes, you could do that. I don't know about stability improvements, but it would definitely be faster.

[Bruno Luong]

But here A is a vector isn't it?

So when we discuss of advantage/disadvantage of

inv(C)

vs

C\eye(size(C))

they are no longer apply to

inv(C)*A 

vs

C\A

For vector C\A win hand down, period.

[Matt J]

We can only assume it is a vector, based on where Armin wrote: "where mu is the centroid of each cluster and x the current feature-vector"

[Christine Tobler]

Based on the formula (x-mu)*cov^(-1)*(x-mu)', you'd expect the output here to be symmetric positive (semi-)definite. If cov is numerically full rank, you can use the following to compute this matrix in a way that will keep it symmetric positive semi-definite.

R = chol(cov);
A = (x-mu)/R;
A = A*A';

If cov is low-rank (positive semi-definite), you could consider doing the same using the ldl function instead of chol.

Answer 1

[Shashank Prasanna]

While inv indeed does involve more operations than \ because it involves inverse as well as multiplication, \ is just one operation. It is safer as you can use it on non-full rank of the system which is determined by QR decomp. I persuade you to look at this literature from Cleve the creator of MATLAB: Linear Equation and of course this link in the documentation.

To support what James said, your example is a very bad one

Not the best but a supporting example is as follows. Using a larger matrix, and making sure the random numbers are the same

rng(5);tic, inv(rand(1000))*rand(1000,1); toc

rng(5);tic, rand(1000)\rand(1000,1); toc

Elapsed time is 0.321691 seconds.
Elapsed time is 0.142591 seconds.

You will have to run this a couple of times. The first run may result in MATLAB performing some optimizations and memory allocations.

[Hunter Colegrove]

Is there a reason inv(A) uses LU decomposition as the algorithm to solve the inverse of the matrix? Is there a computational advantage of some sort?

Answer 2

[James Tursa]

1) Your example is so small as to make the timing irrelevant

2) You burden the \ method with creating another variable, B

3) The \ operation doesn't know ahead of time that B is the identity matrix, so it may not be able to take computational advantage of that fact.

4) Why do you think inv is more computationally intensive than \ ?

[Matt J]

5) Omitting semicolons and displaying to the screen is also adding unnecessary overhead.

6) You have to re-run the code multiple times, preferably inside a function to see a fully optimized execution. There is no way inv(A) on a 5x5 matrix takes 0.006097 seconds on any reasonable machine.

[the cyclist]

Regarding James' point #4, he may have gotten the idea from the documentation of the inv() function:

"In practice, it is seldom necessary to form the explicit inverse of a matrix. A frequent misuse of inv arises when solving the system of linear equations Ax = b. One way to solve this is with x = inv(A)*b. A better way, from both an execution time and numerical accuracy standpoint, is to use the matrix division operator x = A\b."

@Bramen, I think the larger worry with inv() is actually the numerical stability for arrays that are nearly singular. See, for example, this old blog post:

http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/

Answer 3

[Matt J]

It's not that inv(A) is supposed to be more computationally intensive than \. It's that A\b is faster than inv(A)*b, as demonstrated below.

    N=2000;
    A=rand(N);
    b=rand(N,1);

    tic;
     inv(A)*b;
    toc

Elapsed time is 0.509519 seconds.

    tic
     A\b;
    toc

Elapsed time is 0.195016 seconds.

[Daniel]

The least squares estimate in statistics is a very common use of inv(X'*X) and there the size of (X'X) is usually somewhere between 2x2 and 10x10. For those cases inv() appears much faster, as you can see here:

X = normrnd(0,1,100,4);
beta=[1 2 3 4]';
y = X*beta;

tic;
for i=1:10000
inv(X'*X)*X'*y;
end
toc;

tic;
for i=1:10000
(X'*X)\X'*y;
end
toc;
Elapsed time is 0.144074 seconds.
Elapsed time is 0.236353 seconds.

I have to do this operation literally one billion times, so changing the code reduces the time from 3 days to 2 days which is nice... And the constant change recommendation from Matlab is annyoing.

[Daniel]

And btw the QR method is the slowest:

tic;
for i=1:n
[Q,R] = qr(X,0);
inv(R)*inv(R)'*X'*y;
end
toc;

Elapsed time is 0.375314 seconds.

[Matt J]

You shouldn't be comparing with (X'*X)\X'*y nor with an explicit call to qr(). You should be comparing with X\y, which is equivalent (in this case) to the QR method.

    tic;
    for i=1:10000
    inv(X'*X)*X'*y;
    end
    toc;
    tic;
    for i=1:10000
    (X'*X)\X'*y;
    end
    toc;

    tic;
    for i=1:10000
    X\y;
    end
    toc;

Elapsed time is 0.082173 seconds.
Elapsed time is 0.157320 seconds.
Elapsed time is 0.086341 seconds.

That said, when you have a very tall X matrix, then using the normal equations directly can be faster, and also save memory. However, since X'*X has a higher condition number than X, it will also be less stable numerically and amplify noise, as demonstrated in the example below. So you need to be sure you have a well-conditioned X in order for this to be a good compromise.

        d=30;
        N=1e5;

        X=diag(1:d).^2*rand(d,'single');
        noise=randn(d,N);

       z1 = X\noise;
     Error1=std(z1(:)),

       z2=(inv(X.'*X)*X.'*noise);
     Error2=std(z2(:)),

Error1 =

    1.6255

Error2 =

    210.7583

[Daniel]

Oh, thanks for that answer. I wasn't aware of this. It should make the code faster and ease the problem of stability that I happened to have.