how to avoid nan result when calculate sparse matrics division

8 次查看(过去 30 天)
I'm doing a calculation using sparse matrics with the size of 12905*12905. I attached the code here.
WH = ones(P.npts, 1);
A = [L.*WL;sparse(1:P.npts,1:P.npts, WH)];
b = [zeros(P.npts,3);sparse(1:P.npts,1:P.npts, WH)*P.pts];
cpts = (A'*A)\(A'*b);
'WL' is a weighted value. 'P.npts' means the number of points with double type, which is a mumber in the structure 'P'.
P.npts = 12905
My question is: Sparse matric 'A' and 'b'are not a zero matrics, and matric (A'*A) either. But I always get 'NaN' result of cpts matrics. How to avoid getting NaN result when calculating a sparse matric. The structure rank of A is 12905. Could anyone help to solve this problem? Really appreciated!
Thanks

采纳的回答

John D'Errico
John D'Errico 2017-8-31
Sigh. Without seeing your numbers, this is hard to tell. Just telling us the size of the matrix is irrelevant.
I will tell you that you are doing something that is foolish.
Do NOT do this:
cpts = (A'*A)\(A'*b);
That squares the condition number of the problem.
Whereas this does not:
cpts = A\b;
Whoever told you to use the form you are using was steering you to a very poor way of solving the problem.
In the end, what matters is what the rank and condition number of A is in terms of linear algebra. So what does this tell you?
cond(A)
If that is too slow, then just do
condest(A)
It won't be to far off. If the condition number is greater than roughly 1e8, then the form you were using will have a condition number of roughly 1e16, which will then potentially yield NaNs for a result.
So at the very least, report the condition number of A.
  2 个评论
Heather Zhang
Heather Zhang 2017-9-1
Thank you very much, John, for your advise. The size of matric b is 25810*3. The size of sparse matric A is 25810*12905 which is not a square matrics. Could you tell me how to estimate the condition number of A? Rearlly appreciate it. The condition number of (A'*A) is inf. Thanks

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Sparse Matrices 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by