3 ms·

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 po

by 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
    #include <stdio.h> #include <math.h> double interp(double x, double x1, double y1, double x2, double y2) { double node_poly = (x - x1) * (x - x2); double w1 = 1.0 / (x1 - x2); double w2 = 1.0 / (x2 - x1); double y = node_poly * ((w1 * y1) / (x - x1) + (w2 * y2) / (x - x2)); return y; } int main() { for (int p = -300; p >= -330; p -= 1) { double x = pow(10, p); double y = interp(x, 0, 10, 1, 20); printf("x = %.6g, y = %.6g\n", x, y); } return 0; }
    - jacobolus 5y ago
      You implemented the wrong formula. You don’t want to explicitly calculate "node poly". The formula you want is the “second barycentric formula”, (4.2) at the top of page 505 of https://people.maths.ox.ac.uk/trefethen/barycentric.pdf https://people.maths.ox.ac.uk/trefethen/barycentric.pdf Or formula (2.6) halfway down page 2 of https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011_143.pdf https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011... Also cf. Higham (2004) “The numerical stability of barycentric Lagrange interpolation” https://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf https://www.maths.manchester.ac.uk/~higham/narep/narep440.pd... If you want to implement the “first barycentric formula” instead you need to take some additional care to avoid overflow/underflow. As the papers above note, the choice of interpolation nodes makes a big difference. If you have an arbitrary function to approximate on an arbitrary complex domain, take a look at https://arxiv.org/abs/1612.00337 https://arxiv.org/abs/1612.00337
jacobolus 5y ago
> 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. You are correct, it does not. Quite the opposite: the results are excellent when x doesn’t quite equal one of the interpolation nodes, because you end up with a value of S · y / S, where S = «very big but not close to overflowing float» and y = «expected output value». The only potentially problematic case is if S · y overflows but S does not. But floating point numbers have a ton of dynamic range, so unless you are looking at points really, really close to zero, with very large outputs (say, x − x_j = 1e-200, y = 1e+200), there is no problem. And the failure case there is that you end up with Inf instead of 1e200 as an output. It isn’t a problem in practice, and if you are worried about it as a problem in theory, you can take care in choice of interpolation nodes or do some scaling of the output values.

by gjm11 5y ago

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
  #include <stdio.h> #include <math.h> double interp(double x, double x1, double y1, double x2, double y2) { double node_poly = (x - x1) * (x - x2); double w1 = 1.0 / (x1 - x2); double w2 = 1.0 / (x2 - x1); double y = node_poly * ((w1 * y1) / (x - x1) + (w2 * y2) / (x - x2)); return y; } int main() { for (int p = -300; p >= -330; p -= 1) { double x = pow(10, p); double y = interp(x, 0, 10, 1, 20); printf("x = %.6g, y = %.6g\n", x, y); } return 0; }
  - jacobolus 5y ago
    You implemented the wrong formula. You don’t want to explicitly calculate "node poly". The formula you want is the “second barycentric formula”, (4.2) at the top of page 505 of https://people.maths.ox.ac.uk/trefethen/barycentric.pdf https://people.maths.ox.ac.uk/trefethen/barycentric.pdf Or formula (2.6) halfway down page 2 of https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011_143.pdf https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011... Also cf. Higham (2004) “The numerical stability of barycentric Lagrange interpolation” https://www.maths.manchester.ac.uk/~higham/narep/narep440.pdf https://www.maths.manchester.ac.uk/~higham/narep/narep440.pd... If you want to implement the “first barycentric formula” instead you need to take some additional care to avoid overflow/underflow. As the papers above note, the choice of interpolation nodes makes a big difference. If you have an arbitrary function to approximate on an arbitrary complex domain, take a look at https://arxiv.org/abs/1612.00337 https://arxiv.org/abs/1612.00337

jacobolus 5y ago

> 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. You are correct, it does not. Quite the opposite: the results are excellent when x doesn’t quite equal one of the interpolation nodes, because you end up with a value of S · y / S, where S = «very big but not close to overflowing float» and y = «expected output value». The only potentially problematic case is if S · y overflows but S does not. But floating point numbers have a ton of dynamic range, so unless you are looking at points really, really close to zero, with very large outputs (say, x − x_j = 1e-200, y = 1e+200), there is no problem. And the failure case there is that you end up with Inf instead of 1e200 as an output. It isn’t a problem in practice, and if you are worried about it as a problem in theory, you can take care in choice of interpolation nodes or do some scaling of the output values.