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
    Commented Nov 21, 2014 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$ Commented Nov 21, 2014 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$ Commented Nov 21, 2014 at 18:06

1 Answer 1

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
    Commented Nov 21, 2014 at 4:07
  • $\begingroup$ @horchler that's a good point, I've updated my answer a bit. $\endgroup$
    – David
    Commented Nov 21, 2014 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$ Commented Nov 21, 2014 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
    Commented Nov 21, 2014 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$ Commented Nov 21, 2014 at 18:24

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .