What exactly does ε represent in numerical computing

I often see ε (epsilon or flipped 3, whatever you call it) in different contexts and I'm not sure I fully grasp it. Sometimes it's the machine epsilon (eps in MATLAB), sometimes it's just a small positive number in proofs and it's occasionally used in optimization tolerances.

Is machine epsilon specifically the smallest number such that 1 + ε > 1 in floating point? And when people set solver tolerances like 1e-6, is that related to machine epsilon or just an arbitrary choice? Just wanted to hear how people actually think about it 'practically'!

Answer 1

[Stephen23]

Epsilon is just a letter. Like all letters, epsilon can mean lots of different things, just like the letters "alpha" and "gamma" have lots of different meaning depending on the field, the context, the author, the target audience, etc.

"Is machine epsilon specifically the smallest number such that 1 + ε > 1 in floating point?"

That is one common definition, yes. But note that there is no one "floating point", there are many standard (and some esoteric) floating point number systems, all(?) of them have some concept of machine epsilon with different actual values:

https://en.wikipedia.org/wiki/Machine_epsilon#Values_for_standard_hardware_arithmetics

"And when people set solver tolerances like 1e-6, is that related to machine epsilon or just an arbitrary choice?"

Arbitrary choice: they selected a letter and defined it to mean something. You can do that too.

Some letters are used by convention for particular meanings, e.g. pi. But even pi, famously used for the transcendental constant, is also used with plenty of other meanings.

[William]

Fair point, I think I was conflating the formal definition with its various informal uses. Makes more sense now that it's just context dependent like any other symbol. Thanks so much Stephen!

Answer 2

[John D'Errico]

Just to add some more info...

For example, pi(x) is known as the prime counting function to some, the number of primes less than or equal to x.

https://en.wikipedia.org/wiki/Prime-counting_function

Just one of the various alternatives pi can be used to represent. There are too few greek letters for all the things one might want to represent. As a name, eps also has many various meanings, though the most common one in my experience is to represent a small thing. (But then, I've always been primarily a numerical animalist.)

A solver tolerance of 1e-6 is just a relatively small number. One hopes it is small enough that you are sufficiently close to a "solution". One part in a million seems small to many. But depending on whether that tolerance is an absolute or a relative tolerance can be a hugely different thing.

And for SOME (approximately quadratic) problems, one might want to use a tolerance that is at least close to sqrt(eps), where here eps is machine epsilon for the given precision you are working in. As a double, we have

eps('double')

ans = 2.2204e-16

sqrt(eps('double'))

ans = 1.4901e-08

This is at least reasonably close to the arbitrary 1e-6 tolerance many use.

Finally, is eps the smallest number you can add to 1, and see a different result? I think this is a common misperception.

hexnums = ['3ff0000000000000'
        '3ff0000000000001'];

These two numbers are represented in hex, different by one bit, down in the least significant bit in the vicinity of 1.

dubnums = hex2num(hexnums);
num2str(dubnums,55)

ans = 2×54 char array
    '                                                     1'
    '1.0000000000000002220446049250313080847263336181640625'

You should see the first is EXACTLY 1. The second is different by eps('double'), or 2.2204e-16.

diff(dubnums)

ans = 2.2204e-16

So the stride between two consecutive doubles near 1, is eps. A 1 bit difference in the LSB. But is that actually the smallest number we could have added to 1, and still get a different result?

delta = 1;
x = 1;
xhat = x + delta;
d = 1023/1024;
while x ~= xhat
  delta = d*delta;
  xhat = x + delta;
end
delta = delta/d

delta = 1.1103e-16

As you can see, delta here is SMALLER than eps. If I went any smaller, then we had x==x+delta.

x == x+delta

ans = logical
   0

x == x + delta*d

ans = logical
   1

A better way to think of eps is as the relative "stride" between two consecutive floating point numbers.

What happens when we added delta to 1? We get the number 1+1.1103e-16, which is then ROUNDED to the closest double precision number in that vicinity, that is...

num2str(1+1.1103e-16,55)

ans = '1.0000000000000002220446049250313080847263336181640625'

So you should now understand that delta should be eps/2. delta as I created it is 1/2 of the stride between consecutive floating point numbers, and eps is that stride between consecutive floats.

eps('double')/2

ans = 1.1102e-16

For larger numbers of course, that stride is itself larger in absolute value. For example, 1024 is a double precision number. I would expect eps(1024) to be proportionally larger, as it is:

eps(1024)

ans = 2.2737e-13

1024*eps

ans = 2.2737e-13

Of course, in a proof, where epsilon is just used as a small number, in a limit where epsilon goes to zero, etc., there it is just some tiny number. And this has essentially nothing to do with the computationally rigorous definition of eps, in terms of the stride between floating point numbers.

Ask in a comment if I can clarify this anymore.

[Paul]

If eps(x) is the stride between from x to xp, where xp is the next larger floating point number from x, then what is the stride from xn to x, where xn is the next smaller floating point number from x? Is that just eps(x)/2?

For example, with x = 1, we have

[x,xp,xn] = deal('3ff0000000000000','3ff0000000000001','3fefffffffffffff');
x = hex2num(x); xp = hex2num(xp); xn = hex2num(xn);
[isequal(eps(x),xp-x), isequal(eps(x)/2,x-xn)]

ans = 1×2 logical array
   1   1

Does the relationship x - xn == eps(x)/2 hold in general?

[John D'Errico]

@Paul

Good question, but no. Not quite. Consider x0 here.

x0 = 1.25;
hx0 = num2hex(x0)

hx0 = '3ff4000000000000'

Now I'll construct the pair of numbers just 1 LSB on either side of x0 as respectively

hxlo = '3ff3ffffffffffff';
hxhi = '3ff4000000000001';

xlo = hex2num(hxlo); % convert to doubles
xhi = hex2num(hxhi);

x = [xlo,x0,xhi]

x = 1×3

1.2500 1.2500 1.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

format long g
diff(x)

ans = 1×2

2.22044604925031e-16 2.22044604925031e-16

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

As you see, the lower and upper strides were identical. At some rare locations, the stride below a point and the stride above are different, as you showed. In fact, the locations where the stride changes must be rare. In fact, I think we can list the set of locations where the stride below changes from the stride above. If I had to guess, it would be all integer powers of 2. For example...

X = 2.^(-5:5);
hX = num2hex(X)

hX = 11×16 char array
    '3fa0000000000000'
    '3fb0000000000000'
    '3fc0000000000000'
    '3fd0000000000000'
    '3fe0000000000000'
    '3ff0000000000000'
    '4000000000000000'
    '4010000000000000'
    '4020000000000000'
    '4030000000000000'
    '4040000000000000'

Now, construct the pars of doubles just above and below.

hXhi = num2hex(X + eps(X))

hXhi = 11×16 char array
    '3fa0000000000001'
    '3fb0000000000001'
    '3fc0000000000001'
    '3fd0000000000001'
    '3fe0000000000001'
    '3ff0000000000001'
    '4000000000000001'
    '4010000000000001'
    '4020000000000001'
    '4030000000000001'
    '4040000000000001'

Clearly all of those numbers were 1 LSB above. But as I said, this is where the break point appears, the stride just below those powers of 2 is only eps(X)/2.

hXlo = num2hex(X - eps(X)/2)

hXlo = 11×16 char array
    '3f9fffffffffffff'
    '3fafffffffffffff'
    '3fbfffffffffffff'
    '3fcfffffffffffff'
    '3fdfffffffffffff'
    '3fefffffffffffff'
    '3fffffffffffffff'
    '400fffffffffffff'
    '401fffffffffffff'
    '402fffffffffffff'
    '403fffffffffffff'

However, I believe the stride does not change at any other locations besides powers of 2. This is something that deserves better proof than I have offered, but you get what you pay for. ;-) As I showed above, the sub and super strides are the same at almost all points on the number line.

[Paul]

Proof by Wikipedia: Double-precision floating-point format - Wikipedia

"The spacing as a fraction of the numbers in the range from 2^n to 2^(n+1) is 2^(n−52)."

Trying it for a couple of values

isequal(eps(4.5),2^(2-52))

ans = logical
   1

isequal(eps(2^10 + 51),2^(10-52))

ans = logical
   1

isequal(eps(2^-8 + 1e-4),2^(-8-52))

ans = logical
   1

But the the wikipedia page doesn't explain how it works for negative numbers, i.e., does eps(x), with x < 0 take us to the next number away from or towards zero? According to eps it's the former (I think).

Also there must be a separate rule for very small numbers, e.g.

format long
eps(0)

ans = 
    4.940656458412465e-324

format hex
eps(0)

ans = 
   0000000000000001

2^-1023 * 2^-51

ans = 
   0000000000000001

[James Tursa]

In floating point computing, I've always thought of eps(x) as simply the value of the least significant bit of x (not considering the sign), without the need for any discussion of strides to the next or previous representable number or rounding modes.

For the eps(0) case, the 0 bit pattern can be considered just the smallest case of denormalized numbers since the exponent bits are all 0 (the defining characteristic of denormalized numbers), so the value of the least significant bit of 0 is the same as it is for all of the denormalized numbers since denormalized numbers are effectively fixed point.

eps(0) % consider 0 as the smallest denormalized number

ans = 4.9407e-324

eps(ans) % the next value up has same eps

ans = 4.9407e-324

eps(2*ans) % twice the value has the same eps

ans = 4.9407e-324

eps(4*ans) % four times the value has the same eps

ans = 4.9407e-324

etc.

Doesn't matter what value a denormalized number is, the value of the least significant bit is always the same. This pattern will continue up until you multiply by a large enough number to turn the denorm into a normalized number.

eps(2^52*ans)

ans = 4.9407e-324

eps(2^53*ans)

ans = 9.8813e-324

[John D'Errico]

@James Tursa - I like that way of visualizing it. It is easy to think of eps(X) as something that continuously changes, proportional to X itself, because of course we tend to think of floats as real numbers. But the magnitude of that LSB is not a continuous thing, but a very discrete thing. You've made also a good point about eps as applied to denormalized numbers.