Can I find the inverse of a sparse matrix faster?

130 views (last 30 days)
ektor on 15 Apr 2017
Commented: ektor on 27 Apr 2017
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 on 15 Apr 2017
Hi Staf, is this a question about a general matrix NN or just this particular one?
ektor on 27 Apr 2017

John D'Errico on 15 Apr 2017
Edited: John D'Errico on 15 Apr 2017
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.
(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.

David Goodmanson on 16 Apr 2017
Edited: David Goodmanson on 17 Apr 2017
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.