Choosing the right epsilon for comparing doubles

Discussion in 'C++' started by walkietalkie74, Feb 2, 2014.

  1. 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!
     
    walkietalkie74, Feb 2, 2014
    #1
    1. Advertisements

  2. walkietalkie74

    Jax Guest

    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, Feb 2, 2014
    #2
    1. Advertisements

  3. 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).
    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, Feb 2, 2014
    #3
  4. 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...
    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, Feb 2, 2014
    #4
  5. Thanks. That's another option I'd thought about.
    I will google it and see if I can find any libraries.
     
    walkietalkie74, Feb 2, 2014
    #5
  6. walkietalkie74

    Stefan Ram Guest

    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.
     
    Stefan Ram, Feb 2, 2014
    #6
  7. 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.
     
    walkietalkie74, Feb 2, 2014
    #7
  8. walkietalkie74

    Jax Guest

    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, Feb 3, 2014
    #8
  9. walkietalkie74

    red floyd Guest

    red floyd, Feb 4, 2014
    #9
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.