4 ms·

> When the input is one of those points directly, That's against one of the rules of well-behaved floating point programming: Never test for equality.

adrian_b 5y ago
Like any rule, this has exceptions. You should never test for equality numbers which are known only approximately, which is the case for most FP numbers. There are nonetheless cases when certain FP values are known exactly and it is OK to test them for equality. In general the rule does not depend on number representation, it applies equally to floating-point numbers, fixed-point numbers, rational numbers and even large integers in some cases. What counts is whether a number is known exactly or only approximately.
- jlokier 5y ago
  There are exceptions. But the type of barycentric interpolation formula that leads to Inf/Inf at an interpolation point isn't one of those exceptions - it's an example which shows why the rule should be borne in mind. When the input is very close to one of the interpolation points, the calculation becomes ill-conditioned, unsurprisingly as two values are diverging towards Inf, so the result can become incorrect as the input approaches the interpolation point, until it becomes Inf/Inf = Nan at some even closer, but non-zero, distance. Before reaching Nan it can also have an interval where the result is +/-Inf. It depends on the interpolation points. But under some circumstances, it is necessary to recognise when the input is close to an interpolation point, and adjust the formula appropriately.
gjm11 5y ago
It's worth being a bit more explicit about this. The appearance of NaNs doesn't result from breaking the "never test for equality" rule. It happens when the point where you're evaluating the function just happens to exactly match one of the interpolation points. But if you decide to deal with the NaN issue by (1) testing whether your input exactly matches any of the interpolation points, or (2) testing the result for NaN, then you're testing for equality or doing something very similar, and then you may have one of the problems the "never test for equality" is meant to prevent: if bad stuff happens when x exactly equals one of the x_j, then perhaps less-obvious bad stuff happens when x almost exactly equals one of the x_j, and an equality test won't catch that. So, does it? Actually, I think not. Suppose x = xj+h with h very small. Then your interpolant is (h times something of normal size) times (1/h times something of normal size + smaller terms) and although one of those is very small and the other very large, ordinary floating-point arithmetic doesn't care about that. There's no catastrophic cancellation or anything. And floating point has much more dynamic range than precision, so to speak, so until h literally reaches zero you aren't running into trouble. (I guess there might be ways of organizing the calculation that do give you catastrophic cancellation, but if you do it the obvious way then that doesn't happen.) So in this particular case you really could just test for equality or test for NaNs. I think. (For some interesting investigation of other aspects of the numerical stability of this kind of interpolation, see https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011_143.pdf https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011.... But note that the authors of this paper are using the formula for extrapolation or, as they prefer to call it, "numerical analytic continuation".)
- jlokier 5y ago
  > There's no catastrophic cancellation or anything. And floating point has much more dynamic range than precision, so to speak, so until h literally reaches zero you aren't running into trouble. I just tried it to be sure, in C with "double", and organized the calculation the way it's described in the paper you cited, formula 2.1. It failed due to catastrophic rational cancellation close to an interpolation point, depending on the values involved. An equality test was not enough. Approaching the first interpolation point, before it transitioned from the correct result, there were some incorrect Inf results prior to it switching to Nan for an input not equal to any interpolation point. This was trivial 1d interpolation between two points, where the points were (0.0 => 10.0, 1.0 => 20.0).
  - gjm11 5y ago
    Hm, interesting. I also tried some numerical experiments to sanity-check and I got reasonable results right up to 1ulp away from the interpolation points. I tried again using exactly your example -- f(x) = 10+10x, interpolation points at 0 and 1 -- and again everything seemed fine. I wonder what we're doing differently, or whether we happened to test at different points. Here's my code: def interpolate1(f,x,xs,ws): p,s = 1,0 for (xj,wj) in zip(xs,ws): p *= x-xj s += wj*f(xj)/(x-xj) return p*s (it's "interpolate1" because I also implemented the other version of the formula described in the article, which for this purpose behaves very similarly). For your example we want f(x) = 10+10x, xs=[0,1], ws=[-1,1]. As I mentioned before, if you implement the calculation of the polynomial called l(x) in the paper -- the thing my code calls p -- in an unfortunate way, you can get cancellation problems, but for this case I'm not sure I see any way for that to happen. I tried calculating l(x) in a different way that could conceivably cause cancellation issues, aimed at this particular interpolant: def foo(x): p = x*x-x s = -1*10/x + 1*20/(x-1) return p*s but this also gave reasonable results for every value of x I tried other than exactly 0 and exactly 1. What exact values of x are giving you what bad results? What language are you implementing this in, with what data types?*
    - jlokier 5y ago

adrian_b 5y ago

Like any rule, this has exceptions. You should never test for equality numbers which are known only approximately, which is the case for most FP numbers. There are nonetheless cases when certain FP values are known exactly and it is OK to test them for equality. In general the rule does not depend on number representation, it applies equally to floating-point numbers, fixed-point numbers, rational numbers and even large integers in some cases. What counts is whether a number is known exactly or only approximately.

jlokier 5y ago
There are exceptions. But the type of barycentric interpolation formula that leads to Inf/Inf at an interpolation point isn't one of those exceptions - it's an example which shows why the rule should be borne in mind. When the input is very close to one of the interpolation points, the calculation becomes ill-conditioned, unsurprisingly as two values are diverging towards Inf, so the result can become incorrect as the input approaches the interpolation point, until it becomes Inf/Inf = Nan at some even closer, but non-zero, distance. Before reaching Nan it can also have an interval where the result is +/-Inf. It depends on the interpolation points. But under some circumstances, it is necessary to recognise when the input is close to an interpolation point, and adjust the formula appropriately.

gjm11 5y ago

It's worth being a bit more explicit about this. The appearance of NaNs doesn't result from breaking the "never test for equality" rule. It happens when the point where you're evaluating the function just happens to exactly match one of the interpolation points. But if you decide to deal with the NaN issue by (1) testing whether your input exactly matches any of the interpolation points, or (2) testing the result for NaN, then you're testing for equality or doing something very similar, and then you may have one of the problems the "never test for equality" is meant to prevent: if bad stuff happens when x exactly equals one of the x_j, then perhaps less-obvious bad stuff happens when x almost exactly equals one of the x_j, and an equality test won't catch that. So, does it? Actually, I think not. Suppose x = xj+h with h very small. Then your interpolant is (h times something of normal size) times (1/h times something of normal size + smaller terms) and although one of those is very small and the other very large, ordinary floating-point arithmetic doesn't care about that. There's no catastrophic cancellation or anything. And floating point has much more dynamic range than precision, so to speak, so until h literally reaches zero you aren't running into trouble. (I guess there might be ways of organizing the calculation that do give you catastrophic cancellation, but if you do it the obvious way then that doesn't happen.) So in this particular case you really could just test for equality or test for NaNs. I think. (For some interesting investigation of other aspects of the numerical stability of this kind of interpolation, see https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011_143.pdf https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011.... But note that the authors of this paper are using the formula for extrapolation or, as they prefer to call it, "numerical analytic continuation".)

jlokier 5y ago
> There's no catastrophic cancellation or anything. And floating point has much more dynamic range than precision, so to speak, so until h literally reaches zero you aren't running into trouble. I just tried it to be sure, in C with "double", and organized the calculation the way it's described in the paper you cited, formula 2.1. It failed due to catastrophic rational cancellation close to an interpolation point, depending on the values involved. An equality test was not enough. Approaching the first interpolation point, before it transitioned from the correct result, there were some incorrect Inf results prior to it switching to Nan for an input not equal to any interpolation point. This was trivial 1d interpolation between two points, where the points were (0.0 => 10.0, 1.0 => 20.0).
- gjm11 5y ago
  Hm, interesting. I also tried some numerical experiments to sanity-check and I got reasonable results right up to 1ulp away from the interpolation points. I tried again using exactly your example -- f(x) = 10+10x, interpolation points at 0 and 1 -- and again everything seemed fine. I wonder what we're doing differently, or whether we happened to test at different points. Here's my code: def interpolate1(f,x,xs,ws): p,s = 1,0 for (xj,wj) in zip(xs,ws): p *= x-xj s += wj*f(xj)/(x-xj) return p*s (it's "interpolate1" because I also implemented the other version of the formula described in the article, which for this purpose behaves very similarly). For your example we want f(x) = 10+10x, xs=[0,1], ws=[-1,1]. As I mentioned before, if you implement the calculation of the polynomial called l(x) in the paper -- the thing my code calls p -- in an unfortunate way, you can get cancellation problems, but for this case I'm not sure I see any way for that to happen. I tried calculating l(x) in a different way that could conceivably cause cancellation issues, aimed at this particular interpolant: def foo(x): p = x*x-x s = -1*10/x + 1*20/(x-1) return p*s but this also gave reasonable results for every value of x I tried other than exactly 0 and exactly 1. What exact values of x are giving you what bad results? What language are you implementing this in, with what data types?*
  - jlokier 5y ago