Robust computation of the mean of two numbers in floating-point?

Question

Let x, y be two floating-point numbers. What's the right way to compute their mean?

The naive way (x+y)/2 can result in overflows when x and y are too large. I think 0.5 * x + 0.5 * y maybe better, but it involves two multiplications (which maybe is inefficient), and I am not sure if it is good enough. Is there a better way?

Another idea I've been playing with is (y/2)(1 + x/y) if x<=y. But again, I am not sure how to analyze this and prove that it satisfies my requirements.

Moreover, I need a guarantee that the computed mean will be >= min(x,y) and <= max(x,y). As pointed out in Don Hatch's answer, maybe a better way of posing this question is: What's an implementation of the mean of two numbers that always gives the most possible accurate result? That is, if x and y are floating-point numbers, how to compute the floating-point number closest to (x+y)/2? In this case, the computed mean is automatically >= min(x,y) and <= max(x,y). See Don Hatch's answer for details.

Note: My priority is robust accuracy. Efficiency is expendable. However, if there are many robust and accurate algorithms, I would pick the most efficient.

Answer 1

I think Higham's Accuracy and Stability of Numerical Algorithms addresses how one can analyze these types of problems. See Chapter 2, especially exercise 2.8.

In this answer I'd like to point out something that isn't really addressed in Higham's book (it doesn't seem to be very widely known, for that matter). If you are interested in proving properties of simple numerical algorithms such as these, you can use the power of modern SMT solvers (Satisfiability Modulo Theories), such as z3, using a package such as sbv in Haskell. This is somewhat easier than using pencil and paper.

Suppose I'm given that $0\leq x\leq y$, and I'd like to know if $z=(x+y)/2$ satisfies $x\leq z\leq y$. The following Haskell code

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

will let me do this automatically. Here test1 fun is the proposition that $x \leq \mathit{fun}(x,y) \leq y$ for all finite floats $x,y$ with $0\leq x\leq y$.

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

It overflows. Suppose I now take your other formula: $z=x/2+y/2$

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

Doesn't work (due to gradual underflow: $(x/2)\times2 \neq x$, which might be unintuitive due to all arithmetic being base-2).

Now try $z=x + (y-x)/2$:

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

Works! The Q.E.D. is a proof that the test1 property holds for all floats as defined above.

What about the same, but restricted to $x\leq y$ (instead of $0\leq x\leq y$)?

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

Okay, so if $y-x$ overflows, how about $z = x + (y/2-x/2)$?

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

So it seems that among the formulas I've tried here, $x + (y/2 - x/2)$ seems to work (with a proof, too). The SMT solver approach seems to me a much quicker way of answering suspicions about simple floating-point formulas than going through floating-point error analysis with pencil and paper.

Finally, the goal of accuracy and stability is often at odds with the goal of performance. For performance, I don't really see how you can do better than $(x+y)/2$, especially since the compiler will still do the heavy lifting of translating this into machine instructions for you.

P.S. This is all with single-precision IEEE754 floating-point arithmetic. I checked $x \leq x + (y/2-x/2) \leq y$ with double-precision arithmetic (replace SFloat with SDouble), and it works too.

P.P.S. One thing to bear in mind when implementing this in code is that compiler flags like -ffast-math (some forms of such flags are sometimes turned on by default in some common compilers) will not result in IEEE754 arithmetic, which will invalidate the above proofs. If you do use flags that enable, e.g., associative addition optimizations, then there's no point in doing anything other than $(x+y)/2$.

P.P.P.S. I got carried away a little looking only at simple algebraic expressions without conditionals. Don Hatch's formula is strictly better.

Answer 2

First, observe that if you have a method that gives a most accurate answer in all cases, then it will satisfy your required condition. (Note that I say a most accurate answer rather than the most accurate answer, since there may be two winners.) Proof: If, to the contrary, you have an accurate-as-possible answer that does not satisfy the required condition, that means either answer<min(x,y)<=max(x,y) (in which case min(x,y) is a better answer, a contradiction), or min(x,y)<=max(x,y)<answer (in which case max(x,y) is a better answer, a contradiction).

So I think that means your question boils down to finding a most accurate possible answer. Assuming IEEE754 arithmetic throughout, I propose the following:

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

My argument that this gives a most accurate answer is a somewhat tedious case analysis. Here goes:

• Case max(abs(x),abs(y)) >= 1.:

  • Subcase neither x nor y is denormalized: In this case the computed answer x/2.+y/2. manipulates the same mantissas and therefore gives the exact same answer as computation of (x+y)/2 would yield if we assumed extended exponents to prevent overflow. This answer may depend on rounding mode but in any case it's guaranteed by IEEE754 to be a best possible answer (from the fact that the computed x+y is guaranteed to be a best approximation to mathematical x+y, and the division by 2 is exact in this case).

  • Subcase x is denormalized (and so abs(y)>=1):

    answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

  • Subcase y is denormalized (and so abs(x)>=1): analogous.

• Case max(abs(x),abs(y)) < 1.:
  • Subcase the computed x+y is either non-denormalized or denormalized-and-"even": Although the computed x+y may not be exact, it is guaranteed by IEEE754 to be a best possible approximation to the mathematical x+y. In this case the subsequent division by 2 in the expression (x+y)/2. is exact, so the computed answer (x+y)/2. is a best possible approximation to the mathematical (x+y)/2.
  • Subcase the computed x+y is denormalized and "odd": In this case exactly one of x,y must also be denormalized-and-"odd", which means the other of x,y is denormalized with the opposite sign, and so the computed x+y is exactly the mathematical x+y, and so the computed (x+y)/2. is guaranteed by IEEE754 to be a best possible approximation to the mathematical (x+y)/2.

Answer 3

For IEEE-754 binary floating-point formats, exemplified by binary64 (double precision) computation, S. Boldo formally proved that the simple algorithm shown below delivers the correctly rounded average.

Sylvie Boldo, "Formal verification of programs computing the floating-point average." In International Conference on Formal Engineering Methods, pp. 17-32. Springer, Cham, 2015. (draft online)

Since division by two is accurate in binary floating-point arithmetic, provided underflow does not occur, it seems intuitively clear that by choosing one of the two formulas $(x+y)/2$ and $x/2 + y/2$ as appropriate (based on the magnitude of the inputs) once should achieve an accurately-rounded average. Boldo's paper shows that for IEEE-754 binary64 any switchover point $C \in [2^{-967}, 2^{970}]$ will suffice. One might pick $C$ so as to provide the best performance for a particular use case.

This yields the following exemplary ISO-C99 code:

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

In recent follow-up work, S. Boldo and co-authors showed how to achieve best possible results for the IEEE-754 decimal floating-point formats by making use of fused multiply-add (FMA) operations and a well-known precision-doubling building block (TwoSum):

Sylvie Boldo, Florian Faissole, and Vincent Tourneur, "A Formally-Proved Algorithm to Compute the Correct Average of Decimal Floating-Point Numbers." In 25th IEEE Symposium on Computer Arithmetic (ARITH 25), June 2018, pp. 69-75. (draft online)

Answer 4

Although it may not be super-efficient performance-wise, there is a very simple way to (1) make sure none of the numbers are greater than either x or y (no overflows) and (2) keep the floating point as "accurate" as possible (and (3), as an added bonus, even though subtraction is being used, no values will ever be stored as negative numbers.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

In fact, if you really want to go for accuracy, you don't even need to perform the division on the spot; just return the values of min(x, y) and difference which you can use to simplify logically or manipulate later.

Answer 5

Convert to higher precission, add the values there and convert back.

There should be no overflow in the higher precission and if both are in the valid floating point range, the calculated number should be inside too.

And it should be in between them, worst case just half of the larger number if the precission is not sufficient.

Answer 6

Theoretically, x/2 can be computed by subtracting 1 from the mantissa.

However, actually implementing bitwise operations like this is not necessarily straightforward, particularly if you don't know the format of your floating point numbers.

If you can do this, the entire operation is reduced to 3 addition/subtractions, which should be a significant improvement.

Answer 7

I was thinking along the same line as @Roland Heath but can't comment yet, here's my take:

x/2 can be computed by subtracting 1 from the exponent (not the mantissa, subtracting 1 from the mantissa is subtracting 2^(value_of_exponent-length_of_mantissa) from the overall value).

Without restriction of the general case, let's assume x < y. (If x > y, relabel the variables. If x = y, (x+y) / 2 is trivial.)

• Transform (x+y) / 2 into x/2 + y/2, which can be performed by two integer subtractions (by one from the exponent)
  • However there's a lower bound on the exponent depending on your representation. If your exponent already is minimal before subtracting 1, this method will require special case handling. A minimal exponent on x will make x/2 smaller than representable (assuming mantissa is represented with an implicit leading 1).
  • Instead of subtracting 1 from the exponent of x, shift x's mantissa to the right by one (and add the implicit leading 1, if any).
  • Subtract 1 from the exponent of y, if it's not minimal. If it's minimal (y is bigger than x, because of the mantissa), shift the mantissa to the right by one (add implicit leading 1, if any).
  • Shift the new mantissa of x to the right according to the exponent of y.
  • Perform integer addition on the mantissae, unless the mantissa of x has been shifted out completely. If both exponents were minimal, the leading ones will overflow, which is ok, because that overflow is supposed to become an implicit leading one again.
• and one floating point addition.
  • Can't think of any special case here; except for rounding, which also applies to shifting described above.