Choosing the right epsilon for comparing doubles

[walkietalkie74]

Hi all,

Sorry if this question has been already asked in the past (I bet so). I would appreciate if you could point me in the right direction.

I do not generally compare doubles because I am well aware of the round-offerrors that make such a comparison meaningless.

However, I am trying to write my program in a test-driven way, so I'd like to write a tester that checks that the result of a specific calculation is actually zero.

My code is the one below. It's calculating a straight line equation based on a point and a vector (in space, a straight line is the intersection of two planes, hence the getPlane1() and getPlane2()), and is checking whether the obtained equation is correct by calculating arbitrary points that are supposed to be on the same straight line, and checking whether their coordinates solve the plane equations:

bool GeometryTester::runTests() {

Point p(2, -2, 1);
Vector v(2, 4, 2);
StraightLine straightLine = Geometry::getStraightLine(p, v);

bool result = true;

for (double scalar = -10.0; scalar <= 10.0; scalar += 0.11131719) {
Point calcPoint = Geometry::getPoint(p, v, scalar);

double s1 = straightLine.getPlane1().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

double s2 = straightLine.getPlane2().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

result = result &&
MathComparer::equalsZero(s1);

result = result &&
MathComparer::equalsZero(s2);
}

return result;
}

For now, in MathComparer::equalsZero(), I am using a (wrong - I know) direct == comparison between the argument passed in and 0.0, but I'd like toreplace it with the appropriate comparison using an epsilon.

Am I correct in thinking that the epsilon should actually be chosen depending on the current calculation and the accuracy I am expecting from it?

I can see, for example, that the solve() method most of the times returns 0..0 but at times something like 6e-12 (which is absolutely sensible, given the imprecise nature of floating-pint types).

Would it be sensible to pass in an arbitrary epsilon to the equalsZero() function?

I have seen examples of usage of std::numerics<double>::epsilon but I am not sure how it should be used in my case...

Can you help me understand this stuff better?

Thanks!

 

[Jax]

Hi all,

Sorry if this question has been already asked in the past (I bet so). I
would appreciate if you could point me in the right direction.

I do not generally compare doubles because I am well aware of the
round-off errors that make such a comparison meaningless.

However, I am trying to write my program in a test-driven way, so I'd
like to write a tester that checks that the result of a specific
calculation is actually zero.

My code is the one below. It's calculating a straight line equation
based on a point and a vector (in space, a straight line is the
intersection of two planes, hence the getPlane1() and getPlane2()), and
is checking whether the obtained equation is correct by calculating
arbitrary points that are supposed to be on the same straight line, and
checking whether their coordinates solve the plane equations:

bool GeometryTester::runTests() {

Point p(2, -2, 1);
Vector v(2, 4, 2);
StraightLine straightLine = Geometry::getStraightLine(p, v);

bool result = true;

for (double scalar = -10.0; scalar <= 10.0; scalar += 0.11131719) {
Point calcPoint = Geometry::getPoint(p, v, scalar);

double s1 = straightLine.getPlane1().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

double s2 = straightLine.getPlane2().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

result = result &&
MathComparer::equalsZero(s1);

result = result &&
MathComparer::equalsZero(s2);
}

return result;
}

For now, in MathComparer::equalsZero(), I am using a (wrong - I know)
direct == comparison between the argument passed in and 0.0, but I'd
like to replace it with the appropriate comparison using an epsilon.

Am I correct in thinking that the epsilon should actually be chosen
depending on the current calculation and the accuracy I am expecting
from it?

I can see, for example, that the solve() method most of the times
returns 0.0 but at times something like 6e-12 (which is absolutely
sensible, given the imprecise nature of floating-pint types).

Would it be sensible to pass in an arbitrary epsilon to the equalsZero()
function?

I have seen examples of usage of std::numerics<double>::epsilon but I am
not sure how it should be used in my case...

Can you help me understand this stuff better?

Thanks!

Walkietalkie as I'm sure you know there's more than one pair of planes
which can describe a given line..... so are there some special constraints
on the 2 planes you obtain which define the straight line? Just wondering.

Anyways you may want to consider always choosing the error (epsilon) to be
a few orders of magnitude greater than the specified inaccuracy of your
floating-point hardware. Alternatively, and I prefer this, epsilon could
be based on the magnitude of the numbers you are comparing rather than an
absolute constant.

 

[walkietalkie74]

Walkietalkie as I'm sure you know there's more than one pair of planes
which can describe a given line..... so are there some special constraints
on the 2 planes you obtain which define the straight line? Just wondering..

Hi Jax, thanks for your answer.
Yes, I know that there can be infinite pairs of planes which can describe agiven line. This is the actual algorithm used in the method that builds the equations of the planes, given a point and a vector:

StraightLine Geometry::getStraightLine(const Point& p, const Vector& v) {
StraightLine straightLine;

double vx = v.getvx();
double vy = v.getvy();
double vz = v.getvz();

double xp = p.getX();
double yp = p.getY();
double zp = p.getZ();

straightLine.setPlane1(Plane(vy, -vx, 0, vx*yp - vy*xp));
straightLine.setPlane2(Plane(0, vz, -vy, vy*zp - vz*yp));

return straightLine;
}

Do you see any errors?
I used the GeometryTester class in my previous post to test that the equations are correct. Almost all the points that I calculate are solutions for the equations, with the exceptions of a few for which the solve() functions returns something like 0.000000000006....(which should be an approximation of zero - though I agree that it's quite a big error).

Anyways you may want to consider always choosing the error (epsilon) to be
a few orders of magnitude greater than the specified inaccuracy of your
floating-point hardware. Alternatively, and I prefer this, epsilon could
be based on the magnitude of the numbers you are comparing rather than an
absolute constant.

Thanks. I have already discarded the idea of an absolute constant. I am aware of the existence of std::numerics<double>::epsilon and I've found a few examples around showing how to relate it to the magnitude of the numbers involved in the calculations.
Someone suggested to use it the following way:

a and b are "equal" if:
std::abs(a - b) <= std::abs(a) * std::numeric_limits<double>::epsilon

It seems using the same algorithm shown in section 4.2.2 of The Art of Computer Programming (D. Knuth), but it doesn't seem to work.

 

[walkietalkie74]

However, most probably your algorithm has a limited number of floating-
point operations, so it ought to be possible to calculate the maximum
possible error. Alas, finding this can be well more complicated than the
original algorithm so it seems an overkill for a test function (you might
need to write a test for the test, etc). Probably the best approach is to
use some kind of epsilon estimated by gut feeling, but not totally
arbitrary.

Thanks for your answer. The problem with choosing an epsilon dictated by gut feeling is that it would make the tester class necessarily unreliable. I wouldn't be 100% sure that the epsilon would be correct in all circumstances, using any numbers. If a new developer adds a test case which fails, I have no way to figure out whether it's because of the wrong epsilon chosen orthe new test condition itself...

6e-12 seems actually pretty big error unless the numbers in your tests
themselves are pretty large. The zero which is output from your test
function is most probably the result of subtracting two large numbers.
You should actually compare the relative difference of those numbers,
e.g. divide them by each other and if the result differs from 1.0 more
than a few epsilons (std::numeric_limits<double>::epsilon()) then your
algorithm might be either wrong or numerically unstable (the latter is
worse). The exact value of "few" above would be based on gut feeling ;-)

Thanks for the suggestions. Really precious. I am afraid you're right. Thatdifference seems quite a big error to me as well. I might well be doing what you say, because I've analysed the values at the intermediate steps of the algorithm and, for example, I've got something like the following:

38.398935600000016 -
3.3989356000000184 -
35

I will change the algorithm and follow your suggestions. Thanks!

 

[walkietalkie74]

The other alternative would be to replace the floating point either with

an "exact math" number class

Thanks. That's another option I'd thought about.
I will google it and see if I can find any libraries.

 

[Stefan Ram]

a straight line is the intersection of two planes

A floating point number, like say, 10000.1, has an IEEE
hexadecimal representation, here: 0x1.3880ccccccccdp13, for
example. This actually represents a whole range of numbers
from 0x1.3880ccccccccc8p13 to 0x1.3880ccccccccd7p13 or
something like that.

This error of the plane parameters gives an error
(uncertainity) of the planes and of the line parameters. The
points of the line have to be points of the plane within the
limits of those errors. The details of error propagation
(propagation of distributions) are a standard subject of
numerical mathematics, I believe.

And strictly for every triple of a line and two planes,
we do not get a yes/no-answer but a probability that it
is correct.

But for any plane-like object in the physical world, we
cannot measure its parameters up to 15 digits, so the
measurement error usually exceeds this IEEE represenation
error. Possibly, this larger error also has to be taken
into consideration.

 

[walkietalkie74]

This error of the plane parameters gives an error

(uncertainity) of the planes and of the line parameters.

Hi Stefan. Thanks. Yes, I know that the floating point values are inherently imprecise and that any calculation with them can propagate errors. I wouldn't normally compare two doubles for equality. This necessity arises from the fact that I need to write a tester class. My three-dimensional objects will eventually be projected and displayed, so I am not really interested in having a perfect accuracy. I will try to rearrange the algorithm.

 

[Jax]

Hi Jax, thanks for your answer.
Yes, I know that there can be infinite pairs of planes which can
describe a given line. This is the actual algorithm used in the method
that builds the equations of the planes, given a point and a vector:

StraightLine Geometry::getStraightLine(const Point& p, const Vector& v)
{
StraightLine straightLine;

double vx = v.getvx();
double vy = v.getvy();
double vz = v.getvz();

double xp = p.getX();
double yp = p.getY();
double zp = p.getZ();

straightLine.setPlane1(Plane(vy, -vx, 0, vx*yp - vy*xp));
straightLine.setPlane2(Plane(0, vz, -vy, vy*zp - vz*yp));

return straightLine;
}

Do you see any errors?
I used the GeometryTester class in my previous post to test that the
equations are correct. Almost all the points that I calculate are
solutions for the equations, with the exceptions of a few for which the
solve() functions returns something like 0.000000000006....(which should
be an approximation of zero - though I agree that it's quite a big
error).


Thanks. I have already discarded the idea of an absolute constant. I am
aware of the existence of std::numerics<double>::epsilon and I've found
a few examples around showing how to relate it to the magnitude of the
numbers involved in the calculations. Someone suggested to use it the
following way:

a and b are "equal" if:
std::abs(a - b) <= std::abs(a) * std::numeric_limits<double>::epsilon

It seems using the same algorithm shown in section 4.2.2 of The Art of
Computer Programming (D. Knuth), but it doesn't seem to work.

Walkie.... I'm a complete tyro at C++ so I'm not able to comment on your
examples.

Nice to see people like you still use Knuth/TAOCP although it's highly
theoretical and section 4.2.2 only contains a handful of theorems (four I
think) whose results are most probably incorporated into library subroutines.

Why are you using vector geometry for this? If your problem had dynamic
properties (such as aircraft dynamics in 3D) then I could understand. But
something as static as your straight line seems to be crying out for a
Cartesian approach. Just saying!

 

[red floyd]

Comparing floating point values is a tricky business. I have written
an article about the subject:

http://warp.povusers.org/programming/floating_point_caveats.html

And the classic document:

http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html