MATLAB Answers

John Doe
1

Calculating parts of the sparse matrix inverse

Asked by John Doe
on 26 Apr 2013

I have a sparse symmetric matrix H, containing on average 4 non-zero elements per row/col, and a sparse vector M with only two non-zero elements, +1 and -1, in positions t and f.

x is the values of the elements in H(i,j).

n = 3000;
H = sparse(i,j,x,n,n);
M = sparse([t f],1,[1 -1],n,1);

I want to find X = M'*H^-1*M. The lu-factors of H are computed earlier.

I know the following works:

X = M' * (U \ (L \ (P \ M)));

However, I want to go from right to left, in such a way that only two colums of P are "visited" in the calculation of (P \ M), and similar for (L \ ( P \ M)).

Is there a way of doing this "manually"?

Thank you!

  0 Comments

Log in to comment.

1 Answer

Answer by Matt J
on 26 Apr 2013
Edited by Matt J
on 26 Apr 2013

However, I want to go from right to left, in such a way that only two colums of P are "visited" in the calculation of (P \ M), and similar for (L \ ( P \ M)).

Well, the calculation of P\M has little to do with the columns of P. It's the rows that matter. In particular, since P is a permutation matrix, remember that inv(P)=P.'. So, if you wanted to optimize y=P\M, you could instead do

 y=(P(t,:)-P(f,:)).';

As for L\y and U\L\y, I don't think there's much you can do to optimize those stages. L and U have no simplified inverse and also the sparsity of U,L, and y are aready taken into account by mldivide.

x = A(t,f); % Not good for sparse matrices

No idea why you think this is "not good". FIND will not be better, though.

  6 Comments

Thanks again Matt!

I'm doing these calculations thousands of times, the location of the elements of M changes, thus speed is essential. This is what I want to achieve (copy from article):

"When the independent vector is sparse, relative few of the columns of L have to be operated on. This is the so-called "fast-forward" substitution. ... If the solution is required for only a few network nodes, a "fast backward" substitution process can be carried out, operating on only a few of the rows of U."

I'm doing these calculations thousands of times, the location of the elements of M changes, thus speed is essential.

If only M is changing, but not H, then the optimal way to calculate X is to create a matrix MM whose columns are the different M. Then you can process all M in a single call to mldivide

X=sum( MM.*(H\MM) ,1)

This is advantageous because then mldivide doesn't have to refactor H for all the different M. See,

http://www.mathworks.com/matlabcentral/answers/68654-1-why-for-dense-matrices-it-s-not-useful-to-split-b-2-what-are-the-factors-that-determine-the-spe

Also, mldivide is multithreaded, so it is likely that the processing of different M will be split into parallelized tasks.

When the independent vector is sparse, relative few of the columns of L have to be operated on.

It's not true, which could account for the lack of a reference or explanation in the article. I'm assuming, of course, that L has no special structure that you haven't mentioned. To see that it's not true for general L, consider this example

    z=rand(1,n-1);
    M=[1;zeros(n-2,1);-1];
    L=eye(n)+diag(-z,-1);

You can readily verify that the solution to L*y=M, for any given vector z, is

y=[1;cumprod(z(:))]-[zeros(n-1,1);1];

The solution therefore involves all of the columns of L through cumprod(z).

That's not to say that the sparsity of L cannot be exploited, and mldivide certainly does so, but the inversion in general can involve all of L's columns.

I'm sure the articles are correct.

I believe I have missed some special structure that would make this possible. I just noticed the sentence: "The key concept is the "factor path" described in (Sparse Vector Methods, Tinney et. al, 1985)."

I haven't had the time to explore it yet, but hopefully the key to the solution is there.

Log in to comment.


Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Win prizes and improve your MATLAB skills

Play today