Determination of the calculation error in Matlab

Hi guys,

could you please explain to me how to best estimate the size of the calculation error caused by the approximate representation of real numbers by floating-point numbers?

I know that eps(x) increases for larger x-values and tryed to consider this in the calculation of the tolerance frame in the following example. Any ideas how the estimation of the tolerance frame could be improved?

rng('default');
a = rand(5);
b = 100;
ab = a*b;

% Compare without tolerance
ab/b == a % contains zeros
% Compare with tolerance
dum = ab/b;
tf = 2*eps(max([abs(dum(:)); abs(a(:))])); % tolerance frame
abs(ab/b -a) <= tf % no zeros
abs(a-b) >= 2*eps(ab)

[Walter Roberson]

Have you considered using ismembertol() ?

[Zwithouta]

ismembertol is using 1e-12 as default tolerance for double-precision inputs. But the calculation error of Matlab can get much bigger for large numbers. How can I dynamically adjust the tolerance frame so that I never run the risk of getting 'false' as ismembertol output, although the compared matrices A and B are actually identical and only differ due to the calculation error.

Answer 1

[John D'Errico]

It sounds like what you want is an interval arithmetic tool, that would quite accurately predict a bracket around the uncertainty around a function like exp(x), given knowledge of the uncertainty in x, which may be initially roughly eps(x).

My guess is you want to do this for more than just the basic question of what is the uncertainty in exp(x), because you may be doing subsequence computations, and wish to know the final tolerance band around the final result.

There are interval arithmetic tools to be found for MATLAB, though I believe the one I found when I looked about a year ago was not free. I even considered writing one myself, but things can get pretty tricky, especially if the bracket can include zero. I think I decided then that it was not worth my effort for the moderately small set of people asking for it. You can actually have situations where a bracket splits into two disjoint brackets, and a good interval arithmetic code might worry about that circumstance. For example, what is the bracket around x/y, where both x and y lie within the bracket [-delta,delta]? Feel free to choose as small a value for delta as you wish.

A less nasty problem arises from non-monotonic relations. For example, what is the band around cos(x), where x lies in the interval

pi+[-delta,delta]

Simple evaluation of cos at the endpoints of the interval will give the same constant since cos(x) is symmetric around pi. But in fact, the theoretical interval will be more like:

[-1, cos(pi+delta)]

So the computation can be fairly difficult to do, at least if you insist on more accuracy than you can just get from a simple eps(exp(x)). For example, suppose we consider

x = 2.1;

eps(1)
ans =
      2.22044604925031e-16

eps(2.1)
ans =
      4.44089209850063e-16

eps(2.1)/eps(1)
ans =
     2

Really, eps(x) is just telling you the magnitude of the least significant bit in the number x. But, suppose we decide that x lies in the interval [x-eps(x)/2,eps+eps(x)/2]

x = 2.1;
exp(x + eps(x)*[-1 1])
ans =
        8.16616991256765          8.16616991256766

exp(x + eps(x)*[-1 1]) - exp(x)
ans =
      -3.5527136788005e-15      5.32907051820075e-15

Note that I get the exact same result for that as when I reduced the interval width by a factor of 2.

exp(x + eps(x)*[-1 1]/2) - exp(x)
ans =
      -3.5527136788005e-15      5.32907051820075e-15

eps(exp(x))
ans =
      1.77635683940025e-15

Anyway, there really is little you can do here, unless you go to a full fledged interval analysis tool. And even then, if you really insist on true accuracy from these computations, you would need to go to higher precision arithmetic.

For example, it is fairly easy to use a tool like HPF to give at least a simple bracket on the result. Thus, done in 22 digits of precision, with the last two digits hidden as "guard" digits, we could write this:

x = hpf('2.1',[20,2]) + eps(2.1)/2*[-1 1]
x =
HPF array of size: 1  2
|1,1| 2.0999999999999997780
|1,2| 2.1000000000000002220

exp(x)
ans =
HPF array of size: 1  2
|1,1| 8.1661699125676482602
|1,2| 8.1661699125676518867

So the interval around the result was

exp(x) - exp(hpf('2.1',[20 2]))
ans =
HPF array of size: 1  2
|1,1| -1.813253973e-15
|1,2| 1.813253973e-15

What did I say before was the least significant bit on exp(2.1)?

eps(exp(2.1))
ans =
      1.77635683940025e-15

So did I really gain much by this exercise? Not a lot. And remember that while I am quite confidant in the accuracy of the HPF computations (since I know and trust the skills of the author quite well), it is true that HPF will be quite a bit slower than similar computations done with doubles. As well, the simple interval computations I did there were valid for a nice monotonic function like exp(x), but things can get nasty if one is not careful.

You can find HPF on the File Exchange, or I suppose you could have used syms to do similar computations. HPF is free at least.

36534-hpf-a-big-decimal-class

A quick search for "interval arithmetic MATLAB" finds a couple of tools. The first one (b4m) is clearly free, but it looks not to run on a MAC, which is important for me, since I could not test it. The second one (INTLAB) is not free.

http://www.ti3.tu-harburg.de/zemke/b4m/

http://www.ti3.tu-harburg.de/rump/intlab/

Honestly, only you can decide if it is worth the investment of time and possibly money to do this. A year ago, I decided not to bother, and as I said, I know the fellow who wrote HPF. IMHO, a thorough understanding of numerical methods is usually sufficient to be confidant in my results.

Answer 2

[Walter Roberson]

When calculations are done using matrix inversion or equivalent (such as the \ or pinv() operations) then the calculation error can be arbitrarily high as the matrices approach singularity.

When integration is done, then given any particular deterministic numeric integration algorithm, it is possible to construct a function for which the numeric integral will have arbitrarily high error.

Therefore, unless a particular sequence of operations is done, about the best we can offer is to say that if you want to be sure that two arrays will not be falsely determined to be different when algebraically they are the same but differ due to different sequence of operations and numeric round-off, then what you need to do is to declare that all arrays that are the same size are to be considered to be equal (including the arrays that have nan or inf in them.)

Analyzing a particular sequence of operations properly can take grad-level courses and a lot of experience. It is not as simple as just taking the step-by-step errors, because in any sufficiently complicated algorithm, sources of inaccuracy can affect each other. Worst case in one step might involve data that cannot possibly have worst case in another. A good deal of complex proof is involved. Even the people who design proofs of error bounds for the IEEE 754 algorithms can later be found to have made a mistake.

Speaking of IEEE 754: the more fundamental calculations are handed over to the CPU which is assumed to be IEEE 754 compliant. But IEEE 754 does not mandate particular algorithms, only error bounds -- so you can know the maximum permitted error in any of the algorithms. On the other hand, MATLAB does not use IEEE 754 hardware for some of the operations, in order to improve accuracy under some conditions, and for reproducibility between platforms. But still it can help to study IEEE 754 specifications.

Operations on larger matrices might be handed over to high performance multithread libraries, especially BLAS or LAPACK or Intel's MKL . Those can end up having slightly different results on different architectures due to different hardware instructions. AMD does not use MKL. The relevant libraries also get updated between releases, such as to add in additional hardware vectorization.

The error bounds designs for LAPACK are discussed at http://www.netlib.org/lapack/lug/node72.html

Answer 3

[tarokh habibzadeh]

72.43-(65.21+7.22)=1.4210854715202e-14

72.43==(65.21+7.22)

logical 0

why?

[Walter Roberson]

MATLAB does not calculate in decimal, it calculates in finite precision binary. 1/10 is not exactly representable in binary.

[tarokh habibzadeh]

Thanks