2

I am looking for a reason to choose one of the following way to compute geometric mean of a long series of floating point x:

  1. take nth root of each x, then multiply all of them
  2. multiply all of them, then take nth root

I have heard that for floating point numbers, multiplications and divisions lose less information than additions and subtractions do. Therefore I am not considering the sum-exponent trick.

Should I compute geometric mean via 1 or 2, and why?

Update 1, in response to comment:

All x is less than 1 and in double precision. Their order of magnitude ranges between 10^-1 to 10^-6. Please assume the most common method for computing n-th root, since I am using a built-in function of a programming language. Instead of overflow, I am worried about underflow (?) since all values are less than 1. You can assume the length of the series of x to be in order of 10^8

Patrick the Cat
  • 2,138
  • 1
  • 16
  • 33
  • How long is the sequence of floating-point numbers (order of magnitude)? How do you plan to compute the *n*-th root? For method (2), could there be danger of overflowing the floating-point format of your choice? Is performance any consideration? – njuffa Jun 09 '16 at 02:13
  • 3
    Assuming there is no danger of overflowing double-precision floating point format and that accuracy is the most important design criterion, I would suggest using method (2), using a compensated product, then taking the *n*-th root at the end. For the compensated product (which can be computed very efficiently if FMA is available), see: Stef Graillat, "Accurate Floating-Point Product and Exponentiation", *IEEE Transactions on Computers*, Vol. 58, No. 7, July 2009, pp. 994-1000. ([online copy](http://www-pequan.lip6.fr/~graillat/papers/IEEE-TC-prod.pdf)) – njuffa Jun 09 '16 at 04:11
  • 1
    Just a comment: in the exp-sum-log method, it's not really the additions and subtractions that are the main cause of accuracy loss (and in any case there are good algorithms for correctly rounded summation); it's the log calls. Not surprisingly, a function that maps the entire positive floating-point range into the interval (0.0, 710.0) introduces significant loss of information. @njuffa: Any chance of an answer encapsulating your comment? – Mark Dickinson Jun 09 '16 at 08:08
  • Similarly, the nth root operation is going to lose significant information for large `n`; I agree with @njuffa that (2) is going to be significantly more accurate. – Mark Dickinson Jun 09 '16 at 08:11
  • @MarkDickinson I don't want to write answer before the asker had a chance to clarify the items I inquired about. – njuffa Jun 09 '16 at 08:28
  • @njuffa: Seems reasonable; I'd love to see the Graillat link in an answer at some point, though (I wasn't aware of that paper before). Oh, and while I'm here, (0.0, 710.0) should have said (-750.0, 710.0), of course (and I should have specified IEEE 754 binary64 format). It's too early in the morning. – Mark Dickinson Jun 09 '16 at 08:32
  • 1
    The risk of overflow in method 2 can be avoided by accumulating the floating-point exponents separately. –  Jun 09 '16 at 19:54
  • @njuffa All `x` is less than 1 and in double precision. Their order of magnitude ranges between 10^-1 to 10^-6. Please assume the most common method for computing n-th root, since I am using a built-in function of a programming language. Instead of overflow, I am worried about underflow (?) since all values are less than 1. You can assume the length of the series of `x` to be in order of 10^8. – Patrick the Cat Jun 09 '16 at 20:26
  • Is there a reference implementation of geometric mean somewhere, that is very accurate? I have no idea how to judge any of these approaches wrt accuracy, or even basic correctness – jberryman Sep 02 '22 at 21:13

2 Answers2

7

In general, in a sequence of floating-point operations that also involves contracting operations such as square root or cube root, it is advantageous from an accuracy perspective to perform the contracting operations last. For example, sqrt(1.0/x) is more accurate than 1.0/sqrt(x), sqrt(a*b) is more accurate than sqrt(a)*sqrt(b), and cbrt(a*b*c) is more accurate than cbrt(a)*cbrt(b)*cbrt(c).

As a consequence, unless there is a danger of overflowing or underflowing the chosen floating-point format, such as IEEE-754 binary64 (e.g. double in C/C++), in intermediate computation, method [2] should be chosen. Additional aspect relevant to accuracy: if n-th root is computed by exponentiation, such as pow() in C/C++, additional error will be introduced with every computed root, as explained in case of cube root in my answer to this question. Finally, the computation of the n-th root will be slower than a multiplication, so doing only multiplies with a final root computation at the end will also be a superior approach performancewise.

Very accurate results can be achieved with method [2] by using a compensated product (akin to the compensated addition provided by Kahan summation). See the following paper for details:

Stef Graillat, "Accurate Floating-Point Product and Exponentiation", IEEE Transactions on Computers, Vol. 58, No. 7, July 2009, pp. 994-1000 (online)

This compensated product can be computed particularly efficient on systems that provide the FMA (fused multiply-add) operation in hardware. This is the case for all the common modern processor architectures, both CPUs and GPUs. C/C++ provide convenient access to this via the standard math functions fma(), fmaf().

Update: Asker clarified in comment that the risk of underflow is imminent since there are on the order of 108 factors in [10-6, 10-1]. One possible workaround mentioned by @Yves Daoust in a comment is to separate the factors into mantissa and exponent and accumulate them separately. Whether this is practical will depend on the floating-point environment. While C and C++ provide the standard function frexp() for performing this splitting, this function may not be very fast.

Community
  • 1
  • 1
njuffa
  • 23,970
  • 4
  • 78
  • 130
0

Here's what I would do:

  1. The logarithm of the product over a_i is the same as the sum of the logarithms of multiple values of a_i.

  2. The exponent of a logarithm of some value returns the original value.

  3. The n^th root of an exponential function is given by the same exponential function, where the exponent is divided by n.

Putting this all together, the geometric mean of a_i is given by taking the sum of the logarithms of a_i, dividing by n, and taking the exponent.

In short, geom_mean(A) = e^(sum(ln(a_i)) / n), or e^mean(ln(a_i)).

MrAcurite
  • 1