Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Can anybody please kindly explain to me why the members of the output vector produced by the following example are not equal to zero?

mm = 1000;
Clear[i, j, ss, EIG];
ss = ConstantArray[0, {mm, mm}];
%%% "Generating a random sparse matrix (ss)";
For[i = 1, i < (mm) + 1, i++,
  For[j = 1, j < (mm) + 1, j++,
    If[i == j,
     ss[[i,j]] = (RandomReal[{-Mod[{i + j}, 2][[1]],Mod[{i + j}, 3][[1]]}, 1])[[1]];
    If[Abs[i - j] == 1,
     ss[[i,j]] = (RandomReal[{-Mod[{i + j}, 2][[1]],Mod[{i + j}, 3][[1]]}, 1])[[1]]+0.2;
    If[Abs[i - j] == 2,
     ss[[i,j]] = (RandomReal[{-Mod[{i + j}, 2][[1]],Mod[{i + j}, 3][[1]]}, 1])[[1]]-0.1;
    If[Abs[i - j] == 3,
     ss[[i,j]] = (RandomReal[{-Mod[{i + j}, 2][[1]],Mod[{i + j}, 2][[1]]}, 1])[[1]]+0.1;

EIG = Sort[(Select[Cases[(Eigenvalues[ss, -10]), _Real], # > 0 &]),Less];
EIG - Sort[(Select[Cases[(Eigenvalues[ss, -10]), _Real], # > 0 &]),Less]

For instance, one run of this example leaded to this vector:

{8.41341*10^-16, -2.57433*10^-15, -1.9082*10^-16, -8.32667*10^-17}

I am wondering why the output vector members are not zero! You can test this example on your PC. THNAKS.

share|improve this question

closed as too localized by rcollyer, m_goldberg, whuber, Sjoerd C. de Vries, Yves Klett Mar 11 '13 at 17:18

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

They are equal to zero, within numerical tolerances, and that is the key whenever you are using numerical algorithms. – rcollyer Mar 11 '13 at 12:19
Just to elaborate on @rcollyer's comment: You can get rid of the small non-zero terms by using Chop or by using higher precision (see SetPrecision for example). – sebhofer Mar 11 '13 at 12:22
@sebhofer you are absolutely correct, but I would like to add a caution against using Chop indiscriminately. Since your are truncating the numbers, you are effectively adding in rounding errors, if done to early in the calculation, and can seriously impact the results. So, Chop away, but do it at the end. – rcollyer Mar 11 '13 at 12:28
@rcollyer, thank you for your comments, but still I cannot see where the problem is. Whatever we get from Eigenvalues function is stored in EIG vector, but when we take Eigenvalues results from EIG, a slight difference appears. I am wondering what is the source of these small values in the output vector. I have same problem in a big program in which these values are not as small as this example and cannot be resolved by chop. – mak maak Mar 11 '13 at 13:03
@Nasser, if you repeat this example for larger values of mm (like mm=2000), do you get zero again? – mak maak Mar 11 '13 at 13:06

[Too long for a comment...]

If I understand correctly, the claim is this. Taking Eigenvalues[...,-10] multiple times, on the same matrix, gives results that differ by an ULP or so. If this is indeed the claim, then (1) the post does a really good jjob of disguising it, and (2) this probably can happen for a couple of reasons.

One is if the algorithm uses randomized starting values for any iteration. Since I do not see this behavior on Linux machines I'll bet against that. The other has to do with certain subtleties of the OS as to when 64 bit vs 80 bit register arithmetic is used. Some of the remarks here might help to explain this.

--- edit ---

By coincidence I was looking through an old bug report. The gist was that in some Intel MKL library code one can get two different numeric eigen-results that are typically close together. It does seem to depend on byte boundary alignment, which I believe affects whether or how 80-bit floats registers get used.

--- end edit ---

share|improve this answer
I think this is indeed what the post claims. @randomized starting values: It is very interesting for me that you say that, as it confirms a suspicion I had a few years back, regarding eigenvalue computation in Matlab where you could actually see this behavior. I just never followed up on that. – sebhofer Mar 11 '13 at 15:54

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