Find the Smallest Positive Real Number in MATLAB: Machine Epsilon Calculation

[Gengar]

Simple enough - I'm just trying to find the smallest positive real number, ε, such that 1 + ε ≠ 1 in MATLAB (double precision). So the value 'eps' in MATLAB is actually not quite defined this way, and using this program

min = 0;
max = 1;

test = 1;

while test~=(min+max)/2
    test = (min+max)/2;
    
    if 1+test~=1
        max=test;
    end
    
    if 1+test==1
        min=test;
    end
end

yields the value 1.1102x10^(-16), which does satisfy 1+ε ≠1. This value is very close to eps/2, but is not eps/2 (they differ by 2x10^(-32)). And 1+eps/2=1 anyway so that seemingly can't be right. My question is this - is there a published value for MATLAB's machine epsilon, or is my value correct? My value seems quite annoying in that it feels like it should be eps/2, and in my assignment it referenced eps/2 as the machine epsilon. Cheers!

 

[Antiphon]

Matlab uses the native floating point instructions of the processor it's running on. Epsilon will depend on the machine. The MATLAB value might be the largest one of all the supported platforms.

 

[AlephZero]

Simple enough - I'm just trying to find the smallest positive real number, ε, such that 1 + ε ≠ 1 in MATLAB (double precision).

That is not a good definition. You really want the smallest number such that (a) + (εa) ≠ a for any value of a. Actually that's not a complete definition either, since you also need to consider multiplication and division, not just addition and subtraction.

For a computer CPU with IEEE 64 bit floating point arithmetic, the value of ε is 2^-52 = approximately 2.2044e-16. Your code is probably calculating the floating point number that is closest to, but not equal to, ε/2 (which is where your difference of 10^-32 is coming from).

A simple algorithm (in C) is

int diff(double a, double b)
{
  return (a != b);
}

double eps = 1.0;
while (diff(1.0 + eps, 1.0))
   eps /= 2;
eps *= 2;

Using a separate function to compare the numbers defeats optimizng compilers that figure out (wrongly!) that testing 1+eps against 1 is the same as testing whether eps is zero.

 

[jhae2.718]

In MATLAB, eps is a function that returns machine epsilon. For me it is 2.2204e-16.

If you run edit eps you get:

%EPS  Spacing of floating point numbers.
%   D = EPS(X), is the positive distance from ABS(X) to the next larger in
%   magnitude floating point number of the same precision as X.
%   X may be either double precision or single precision.
%   For all X, EPS(X) is equal to EPS(ABS(X)).
%
%   EPS, with no arguments, is the distance from 1.0 to the next larger double
%   precision number, that is EPS with no arguments returns 2^(-52).
%
%   EPS('double') is the same as EPS, or EPS(1.0).
%   EPS('single') is the same as EPS(single(1.0)), or single(2^-23).
%
%   Except for numbers whose absolute value is smaller than REALMIN,
%   if 2^E <= ABS(X) < 2^(E+1), then
%      EPS(X) returns 2^(E-23) if ISA(X,'single')
%      EPS(X) returns 2^(E-52) if ISA(X,'double')
%
%   For all X of class double such that ABS(X) <= REALMIN, EPS(X)
%   returns 2^(-1074).   Similarly, for all X of class single such that
%   ABS(X) <= REALMIN('single'), EPS(X) returns 2^(-149).
%
%   Replace expressions of the form
%      if Y < EPS * ABS(X)
%   with
%      if Y < EPS(X)
%
%   Example return values from calling EPS with various inputs are
%   presented in the table below:
%
%         Expression                   Return Value
%        ===========================================
%         eps(1/2)                     2^(-53)
%         eps(1)                       2^(-52)
%         eps(2)                       2^(-51)
%         eps(realmax)                 2^971
%         eps(0)                       2^(-1074)
%         eps(realmin/2)               2^(-1074)
%         eps(realmin/16)              2^(-1074)
%         eps(Inf)                     NaN
%         eps(NaN)                     NaN
%        -------------------------------------------
%         eps(single(1/2))             2^(-24)
%         eps(single(1))               2^(-23)
%         eps(single(2))               2^(-22)
%         eps(realmax('single'))       2^104
%         eps(single(0))               2^(-149)
%         eps(realmin('single')/2)    2^(-149)
%         eps(realmin('single')/16)   2^(-149)
%         eps(single(Inf))             single(NaN)
%         eps(single(NaN))             single(NaN)
%
%   See also REALMAX, REALMIN.

%   Copyright 1984-2006 The MathWorks, Inc.
%   $Revision: 5.9.4.8 $  $Date: 2006/04/25 07:18:35 $
%   Built-in function.

 

[Gengar]

Matlab uses the native floating point instructions of the processor it's running on. Epsilon will depend on the machine. The MATLAB value might be the largest one of all the supported platforms.

Agreed, but I'm only interested in macheps for my computer.

That is not a good definition. You really want the smallest number such that (a) + (εa) ≠ a for any value of a. Actually that's not a complete definition either, since you also need to consider multiplication and division, not just addition and subtraction.

For a computer CPU with IEEE 64 bit floating point arithmetic, the value of ε is 2^-52 = approximately 2.2044e-16. Your code is probably calculating the floating point number that is closest to, but not equal to, ε/2 (which is where your difference of 10^-32 is coming from).

Would it be wrong to factorise a + εa = (1+ε)a? In which case the problem is unchanged - and for this task I don't need to take multiplication/division into account. If my factorisation is incorrect, then what would be the smallest value, ε, be such that (a) + (εa) ≠ a for any a? It is not eps=2.2044e-16, as my value still satisfies the above condition. Interestingly, eps/2 does satisfy the condition for a > 2 (integer-wise) but not for a=1 or a=2.

 

[AlephZero]

Would it be wrong to factorise a + εa = (1+ε)a?

That is wrong, unless you can prove they are equal for the computer arithmetic you are using. Unlike mathematics, in general computer arithmetic is not commutative or associative.

For example if the value of a is close to zero, εa may be too small to represent as a number different from 0, even if ε was a number as "big" as 0.1. In that case computing a + εa would give a, but computing (1+ε)a would give a different result.