# 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?

• If this is the wrong StackExchange network for the question, please direct me to where this question should be posted. – Noctis Skytower Jun 14 '16 at 13:28
• What have you tried? Where did you get stuck? We do not want to just do your (home-)work for you; we want you to gain understanding. However, as it is we do not know what your underlying problem is, so we can not begin to help. See here for a relevant discussion. If you are uncertain how to improve your question, why not ask around in Computer Science Chat? You may also want to check out our reference questions. – adrianN Jun 14 '16 at 14:06
• I very much doubt that they returned the same value very often. How exactly did you determine that the values were the same? – gnasher729 Nov 6 '19 at 21:04

"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:

1. Multiplying a finite number by zero returns a zero.
2. Multiplying a finite number by one returns that number.
3. Adding zero to any finite number returns that number.
4. 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.

• I agree, (1-x)*y + x*z is the right way to do it. Although I recently noticed that this lacks a desirable property: It's not necessarily true that interpolate_b(x, y, z) == interpolate_b(1-x, z, y). since if x is epsilon/2, 1-x == 1.0. The solution I found was (1-x) * y + (1-(1-x)) * z, although I am paranoid that my compiler will second-guess me. – Ben Mar 22 '18 at 14:39
• Another proble is that interpolate_b is non-monotonic – Eric Aug 19 '18 at 16:29

According to the document "P0811R2: Well-behaved interpolation for numbers and pointers" (the link @Eric provided), it depends:

1. a+t*(b-a) does not in general reproduce b when t==1, and can overflow if a and b have the largest exponent and opposite signs.
2. t*b+(1-t)*a is not monotonic in general (unless the product ab≤0).

The author continues by providing a short list of desirable properties, before suggesting an implementation that fulfills them all.

It is a worthy read. I suggest creating the method one place accessible to all your projects (e.g. your company's single shared math-library wrapper), and reuse it whenever linear interpolation is needed.