Automatic sparsification algorithm and zero-threshold
8 次查看(过去 30 天)
显示 更早的评论
Matlab uses an automatic sparsification algorithm to convert an array from double to sparse and also to propagate sparsity during arithmetic operations. There is a sparsity threshold value to decide if an element should be regarded as 0. Unfortunately this threshold seems that can't be changed from the user. Also it is really small: ~1e-324 (even smaller than realmin=~1e-308).
e.g,
A=bucky; % non-zeros are 1
subplot(131), spy(A), title('A')
subplot(132), spy(A*1e-323),title('A*1e-323')
subplot(133), spy(A*1e-324),title('A*1e-324')
This has an effect on the propagation of sparsity after arithmetic operations on the initial sparse array. Although some resulted arrays are practically sparse with many values smaller that EPS, these values are stored as non-zeros and the memory requirements explode.
So I guess my question is the following: Is there any way to change that tinny threshold?
Thank you for your time, Petros
9 个评论
James Tursa
2014-6-20
@Petros: realmin = 2.2251e-308 is the smallest full precision double number in IEEE, but it is not the smallest number in IEEE double. Below realmin there is a range of denormalized numbers available with less precision (all the way down to only 1 bit of precision). The smallest denormalized number is 4.9407e-324. All of the normalized and denormalized numbers are non-zero, so they propagate in the sparse matrix calculations. I don't know how to change the behavior from strict non-zero to a tolerance. There is an FEX submission that can clean sparse matrices of small values here:
回答(1 个)
Matt J
2014-6-24
编辑:Matt J
2014-6-24
About the only generic thing I can think of is to do the mldivide operation A\B column-by-column and apply the tolerance to each result. Then, at the very end, recompose the matrix:
[ma,na]=size(A);
[mb,nb]=size(B);
I=cell(1,nb); J=I, S=I;
for j=1:nb
tmp=A\B(:,j);
idx=find(abs(tmp)>tol).';
I{i}=idx;
J{i}(1:length(idx))=j;
S{i}=tmp(idx);
end
result=sparse([I{:}], [J{:}], [S{:}], na,nb);
There are potentially more efficient ways to do the intermediate mldivides A\B(:,j) by pre-decomposing A. Depending on A's structure, for example, an LU decomposition is used, as described here, and that can be done once prior to the loop.
Again, though, if you have enough foreknowledge about the solution to make such truncation both useful and safe, you should probably describe to us what it is about the problem structure that makes that foreknowledge possible. It's the kind of foreknowledge that often leads to good customized solutions.
0 个评论
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!