# Weed out numerical artifacts from matrix inversion

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?

• What is "large" in terms of the dimensions of your sparse matrices? Commented Nov 21, 2014 at 4:36
• 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? Commented Nov 21, 2014 at 10:36
• 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. Commented Nov 21, 2014 at 18:06

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.
• 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. Commented Nov 21, 2014 at 4:07
• @horchler: good point. It seems you say one should ignore values less than $10^{-16}$ times the mode of the large values. Commented Nov 21, 2014 at 4:25
• @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. Commented Nov 21, 2014 at 4:35