Calculating parts of the sparse matrix inverse

2 次查看（过去 30 天）
John Doe
John Doe 2013-4-26
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!

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2013-4-26
编辑：Matt J 2013-4-26
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 个评论
显示 4更早的评论
Matt J
Matt J 2013-4-27
编辑：Matt J 2013-4-27
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
John Doe 2013-4-28
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by