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-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!

[Richard Damon]

On 2/2/14, 11:35 AM, walkietalkie74 wrote:
> 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.

...

>
> 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!
>

Yes, you are correct that "How close to zero is ok" is a question that
is VERY much dependent on the problem at hand. There are ways to analyze
the computation to see how much error would be "expected", to help you
set the bounds.

The other alternative would be to replace the floating point either with
an "exact math" number class (assuming the answers are expressible with
that or a number type that keeps track of possible error (either as an
interval of possible values, or as a estimated value and error bound).

[Jax]

walkietalkie74 <a.laf...@gmail.com> wrote in
news:60d97607-ea5c-4001...@googlegroups.com:

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.

--
Jax :)

[Paavo Helde]

walkietalkie74 <a.laf...@gmail.com> wrote in
news:60d97607-ea5c-4001...@googlegroups.com:

> 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.

[...]

> 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?

Yes. As a trivial example, adding-subtracting an unbounded amount of
floating-point numbers can result in arbitrily large shifts from the
correct value.

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.

>
> 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).

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 ;-)

hth
Paavo

[walkietalkie74]

On Sunday, 2 February 2014 21:42:39 UTC, Jax wrote:

> 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 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).

> 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]

On Sunday, 2 February 2014 21:42:39 UTC, Paavo Helde wrote:

> 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 or the 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. That difference 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]

On Sunday, 2 February 2014 21:15:49 UTC, Richard Damon wrote:

> 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.

[walkietalkie74]

On Sunday, 2 February 2014 23:26:54 UTC, Stefan Ram wrote:

> 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.

[Richard Damon]

One question, is how good does the results NEED to be? You really only
need to check that it is that accurate.

Another way to judge, if this error term corresponds to a distance, then
one reasonable way to set allowable error is reasonable epsilon based on
the magnitude of the input points.

[Jax]

walkietalkie74 <a.laf...@gmail.com> wrote in
news:436a61b6-c1ed-4516...@googlegroups.com:

> On Sunday, 2 February 2014 21:42:39 UTC, Jax wrote:
>
>> 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

>> constraint s on the 2 planes you obtain which define the straight line?

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!

--
Jax :)

[Juha Nieminen]

walkietalkie74 <a.laf...@gmail.com> wrote:
> 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.

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

--- news://freenews.netfront.net/ - complaints: ne...@netfront.net ---

[red floyd]

On 2/4/2014 12:22 AM, Juha Nieminen wrote:
> walkietalkie74 <a.laf...@gmail.com> wrote:
>> 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.
>
> 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