Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

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

John Doe

Products

No products are associated with this question.

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

John Doe on 27 Apr 2013

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."

Matt J on 27 Apr 2013

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.

John Doe on 28 Apr 2013

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.

Matt J

Contact us