Can I find the inverse of a sparse matrix faster?

Dear all,

I have the following code

 T=2000;
 NN = speye(T) + sparse(2:T,1:(T-1),2*ones(1,T-1),T,T);
        u = NN\eye(T);
zz = u'.*u;

Is there a faster way to obtain zz? I suspect that the creation of u slows down the code.

Thank you

[David Goodmanson]

Hi Staf, is this a question about a general matrix NN or just this particular one?

[ektor]

Thanks David, It is just about this particular one.

Answer 1

[John D'Errico]

I'm confused. I already showed in your last question that creating u speeds the code up, as otherwise you would need to form u twice. So why do you think that slows things down?

When you perform the operation (NN\eye(T)), MATLAB stores it somewhere in memory. Clearly, it has to stuff those numbers somewhere.

In your first post, you had this operation:

(NN\eye(T))' .* (NN\eye(T))

So MATLAB had to generate that solution TWICE, storing the result in a pair of temporary arrays! The numbers are still stored, but just stored in a place that will be dumped in the bit bucket as soon as they are used.

In this form,

u = NN\eye(T);
zz = u'.*u;

u is formed ONCE. It is stored in memory as u, then used to create zz.

The same operations are again done, but you save the second call to backslash.

timeit(@() NN\eye(T))
ans =
     0.048416

u = NN\eye(T);
timeit(@() u'.*u)
ans =
     0.032393

Compared to the original form you were using:

timeit(@() (NN\eye(T))' .* (NN\eye(T)))
ans =
      0.12069

As you can see, we can break that total time down as:

0.048416*2 + 0.032419
ans =
      0.12925

Which is pretty close to what we saw before.

The only issue is if there is a faster way to solve this? There might be. At least it is worth considering. But before you go any further, STOP! Look at what you are trying to do. Think seriously about what I will say below.

NN is a large sparse matrix. But the computation you are doing makes no numerical sense at all.

rank(full(NN))
ans =
        1999

So NN has NO inverse. That inverse matrix you were forming as NN\eye(T) is numerical garbage.

svd(full(NN))
ans =

              3
              3
              3
              3
              3
              3
              3
         2.9999
         2.9999
         2.9999
         2.9999

...

         1.0003
         1.0002
         1.0002
         1.0002
         1.0001
         1.0001
         1.0001
              1
              1
              1
              1
              0

Do you see that ZERO at the very end? That matrix is flat out singular. It is rank deficient. It has NO inverse. You cannot invert a singular matrix. Period. If you try to do so, you will get numerical garbage.

cond(full(NN))
ans =
   Inf

Do you not believe me that the result is garbage?

norm(u)
ans =
   Inf

nnz(isinf(u))
ans =
      476776

numel(u)
ans =
     4000000

So more than 10% of the elements of u are +/- inf.

You need to resolve the issue of why it is singular, before you bother to try to speed up code that yields something meaningless.

Answer 2

[David Goodmanson]

Hi Staf, Before attempting something as large as T = 2000 it's never a bad idea to try the same thing for T just a bit smaller, with the advantage that you can use full instead of sparse matrices. Say, T = 10 (also T = 11 to make sure that there is no significant difference between even and odd).

The NN matrix has ones on the main diagonal and twos on the first lower subdiagonal. If you look at its inverse you will start to see some big numbers. Taking

u = inv(NN);
u(T,1)          % lower left corner element

it is not hard to find out that the lower left element equals 2^(T-1) up to a sign. Although technically NN is nonsingular at T = 2000, John is perfectly right that in the real world it is numerically singular and uninvertible. Already by T = 200 the value of u(T,1) is around 8e59 and you can't do much useful with that. Also, the smallest svd value of NN is falling like const x (1/2)^T, so way before T = 2000 Matlab will call it zero.

Is it true that you want element-by-element multiplication of u' and u as opposed to normal matrix multiplication? If it's the former, then since NN is [1] lower triangular with [2] ones on the main diagonal, its inverse u has the same properties as well. The conjugate matrix, u', is upper triangular with property [2]. If you take the element-by-element product of u and u' you always end up with

zz = the identity matrix, eye(T).

That's verified experimentally.