Re: left-right (up-down) parsing may lead to over/underflow of prod

[Shai Ayal]

On Tue, Apr 22, 2008 at 3:00 PM, Rolf Fabian <address@hidden> wrote:
>
>  I consider the following 'feature' as a
>  major drawback of Octave because
>  it infringes the expected commutativity
>  and associativity while building products.
>
>  Parsing products from left to right (for vectors)
>  or up-to-down (columnwise product of matrices)
>  may lead to avoidable over /underflow if
>  some involved factors are extraordinarily
>  large and/or small
>
>  :> computer, version
>  i686-pc-msdosmsvc
>  ans = 3.0.0
>
>  EXAMPLE 1 (prod of a vector)
>  :> x1 = [1e150, 1e160, 1e150, 1e160, 1e-310, 1e-310]
>  It is obvious that the correct product of x is equal to 1,
>  but Octave either overflows if the product is
>  evaluated from left-to-right:
>  :> prod (x)
>  ans = Inf
>
>  or underflows (without any warning)
>  if the product is evaluated from right-to-left
>  :> prod  (fliplr (x))
>  ans = 0
>
>  Thus simply rearranging the order of factors changes
>  Octave's result of the product from 'Inf' to '0', which
>  is unexpected because  generally a product should
>  not depend on the ordering of its factors.
>
>  Expected correct answer is returned by
>  :> prod_via_sum_of_logs = exp (sum (log (x)))
>  ans =  1
>
>  Similar artefacts occur if input to prod is a matrix like:
>  :> y = [1e300, 1e-200; 1e300, 1e-200; 1e-200, 1e300; 1e-200, 1e300]
>
>  :> prod(y)
>  ans =
>    Inf     0
>
>  :> prod (flipud (y))
>  ans =
>      0   Inf
>
>  The correct result can again be obtained via sums of logarithms:
>  :> prod_via_sum_of_logs  = exp (sum (log (y)))
>  ans =
>    1.0000e+200   1.0000e+200
>
>  Of course, some spurious complex noise might appear if some factors
>  are negative while using sums of logarithms, but this can be adjusted
>  easily.
>
>  In order to avoid similar artefacts as outlined above I propose to implement
>  the 'prod' function via (column) sums of logarithms of factors.
>

I see 2 argument against this (but I'm not sure how strong they really are)
1. speed -- you need to calc the logs of the N numbers you are
planning to multiply. I am not sure exactly how slow it is, but surely
much slower than a simple product
2. While in the test cases you give the method you suggest is clearly
superior, I'm quite sure there are some examples in which you loose
precision because of summing logs (I'm too lazy to think this out to
the end, but surely someone can think for me ?).

This looks like one of the many cases where numerical calculations
fail because of limited precision, and you just have to know your data
in order to make it work -- a lot of us call it "job security" :)

Shai