2
$\begingroup$

I am working with the inverses to a set of large sparse matrices (in Matlab). A key indicator for my application is the number of non-zero entries in each row, and I recently discovered that I was getting false results due to numerical problems-- i.e. a value on the order of $10^{-30}$ is still non-zero.

To solve this problem, I decided to screen out numerical artifacts by removing entries in the inverted matrix below a threshold value. I figured I'd detect the threshold by counting up the number of matrix entries in each decade. Where T is the inverse (sparse) matrix:

%% matlab
A = abs(T);    
f = full(floor(log10(A(find(A)))));
u = unique(f);
bins = histc(f,u);

Now, when I look at bins I see that my matrices all indeed have bimodal distributions, with the major peak varying in the vicinity of -5 (i.e. $10^{-5}$) and another peak around -35 (i.e. $10^{-35}$).

But, picking a threshold is problematic. With one matrix there's about a 12-decade gap between the peaks (i.e. no entries at all where $10^{-28} < |a| < 10^{-16}$); whereas with another there is no decade without any entries (though there's still two humps in the distribution).

How should I algorithmically choose a threshold? Should I choose the interior minimum? pick a flat value? is this a good approach? Frankly, even picking an interior minimum is annoying and problem-prone, because of potentially noisy data.

Any feedback?

$\endgroup$
3
  • $\begingroup$ What is "large" in terms of the dimensions of your sparse matrices? $\endgroup$ – horchler Nov 21 '14 at 4:36
  • $\begingroup$ How do you invert large and sparse matrix? If you manage to compute the inverse, the matrix is not actually very large. What do you want to do with that approximate inverse? Do you know about sparse approximate inverses? $\endgroup$ – Algebraic Pavel Nov 21 '14 at 10:36
  • $\begingroup$ Ranging from 1000-5000 elements square; The non-inverse (obverse?) matrices are highly sparse, having around 0.5% nonzero elements. The inverse matrices are somewhat sparse, being permutable to block-triangular form and having maybe 30-50% occupancy. Full disclosure: I presently have only 3 specimens to work with. $\endgroup$ – Brandon Kuczenski Nov 21 '14 at 18:06
1
$\begingroup$

You can test this in Matlab using the eps command. eps is the difference between one and the next largest number Matlab can store. For most systems, eps$=2^{-52}\approx2.22\times10^{-16}$. This suggest that anything below about $10^{-15}$ should be treated as zero. Per horchler's suggestion, if you know an approximate magnitude, N, of your non-zero solutions, you can use eps(N) to get a more relevant number, you could then set, for example, anything less than 10*eps(N) to zero.

If you have lots of small roots, you will have to think about using a different method, or using some fancy techniques to avoid bad floating point errors.

$\endgroup$
5
  • $\begingroup$ This assumes that the size of the large non-zero values of the OPs matrix is around one, but you don't know this. They could be much larger or much smaller. If anything, the function form of eps (a.k.a. unit in last place (ULP)) should be used, but this too may be too naïve. $\endgroup$ – horchler Nov 21 '14 at 4:07
  • $\begingroup$ @horchler that's a good point, I've updated my answer a bit. $\endgroup$ – David Nov 21 '14 at 4:13
  • 1
    $\begingroup$ @horchler: good point. It seems you say one should ignore values less than $10^{-16}$ times the mode of the large values. $\endgroup$ – Ross Millikan Nov 21 '14 at 4:25
  • 1
    $\begingroup$ @RossMillikan: it may or may not work depending on the original matrix and how the inverse is to be used. I'd want to check how good the corrected inverse is. Also, recall that eps relates to adding one number to another, not multiplication. I'd try to avoid inverting the matrix in the first place. $\endgroup$ – horchler Nov 21 '14 at 4:35
  • $\begingroup$ Thanks for the suggestion. I am not sure if I can do that (avoid inversion). What I would really like to do is find cycles in a large directed graph. I am reading about Tarjans algorithm now. $\endgroup$ – Brandon Kuczenski Nov 21 '14 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.