4 ms·

This depends how you define “gone wrong”. In evaluating rational functions (a useful tool for approximating all sorts of other functions), one efficient and wel

This depends how you define “gone wrong”. In evaluating rational functions (a useful tool for approximating all sorts of other functions), one efficient and well behaved algorithm is the “barycentric formula”, based on interpolating function values at various specific input points. When the input is one of those points directly, this formula results in Inf / Inf = NaN. Evaluation code needs to check for this case, and then replace the output by the appropriate value for the input. ±Inf is a very natural output for many kinds of numerical codes.

jlokier 5y ago
> 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).
twic 5y ago
What is the rationale for liking rational functions? Years ago i had some interpolation to do, and used the barycentric rational from Boost [1]. The resulting curve was all over the place. A natural cubic spline gave a much more subjectively sensible-looking fit. [1] https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/barycentric.html https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/mat...
- jacobolus 5y ago
  Without knowing a lot more about the context and what precisely you did, it’s impossible to comment about your negative experience. > What is the rationale for liking rational functions? They are cheap to evaluate, (relatively) easy to prove things about symbolically, and will efficiently approximate any reasonably smooth function.

jlokier 5y ago

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

twic 5y ago

What is the rationale for liking rational functions? Years ago i had some interpolation to do, and used the barycentric rational from Boost [1]. The resulting curve was all over the place. A natural cubic spline gave a much more subjectively sensible-looking fit. [1] https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/barycentric.html https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/mat...

jacobolus 5y ago
Without knowing a lot more about the context and what precisely you did, it’s impossible to comment about your negative experience. > What is the rationale for liking rational functions? They are cheap to evaluate, (relatively) easy to prove things about symbolically, and will efficiently approximate any reasonably smooth function.