Numerical stability of linear interpolation

Is one of these two functions more numerically stable than the other?

def interpolate_a(x, y, z):
return y + (z - y) * x

def interpolate_b(x, y, z):
return y (1 - x) + z x


If so, which one is better and why? After testing these functions with several million inputs, they have always returned the same value. My concern is that with certain inputs, they may return different values, and one will be more correct that the other. Is this possible?

"Numerical stability" is a much vaguer term than most people realise. We typically use it when referring to an approximation method, such as some kind of linear analysis, or numeric quadrature, or solving (possibly partial) differential equations. It refers to the property that given some appropriate assumptions (e.g. the inputs are reasonable), the algorithms converges to the ideal solution reasonably, even in the presence of floating-point error.

That's actually not what you're referring to here. However, there is a good answer to the question that you are trying to ask. If x is between 0 and 1 inclusive (and this is an important proviso; it's called "interpolate" but could be used to extrapolate), then interpolate_b is better.

In any sensible system of floating point (e.g. IEEE-754), you can assume several things:

• Multiplying a finite number by zero returns a zero.

• Multiplying a finite number by one returns that number.

• Adding zero to any finite number returns that number.

• For any finite number x, x - x returns a zero.

Standard provisos: Inf is not a finite number. NaN is not a number at all. Negative zero is a thing.

The negative zero thing is occasionally relevant, so when I say "returns a zero", or "returns that number" when that number happens to be zero, do be aware that you may get a zero of the opposite sign that you might expect under some circumstances. The exact rules are complicated, and we won't go into them here.

Still, it's more correct to say, for example, that if x is a finite number which is not zero, then x+0 and x+(-0) are both identical to x. If x is zero, then x+0 and x+(-0) may not be identical, and may not be identical to x, but they will differ from x by at most a sign bit.

First, consider the case where x is zero (and neither y nor z are negative zero). Then both interpolate_a and interpolate_b return exactly y. So far, so good.

Now consider the case where x is exactly 1. Then interpolate_a returns y + (z - y), because the multiply by x is a no-op. This may not be exactly equal to z due to round-off error. (Exercise: Find two single-precision IEEE-754 numbers which have this property in the default rounding mode.)

However, if you follow through the logic above, you can see that if x is exactly 1 (and neither y nor z are negative zero) then interpolate_b returns exactly z.

It's almost always considered a desirable property of any interpolate function that they return exact values at the endpoints. It's this property that makes interpolate_b a better choice.

StackExchange Computer Science Q#59625, answer score: 6