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

[Jaroslav Hajek]

On Tue, Apr 22, 2008 at 2:33 PM, Rolf Fabian <address@hidden> wrote:
>
>  -------- Original-Nachricht --------
>  > Datum: Tue, 22 Apr 2008 14:12:12 +0200
>  > Von: David Bateman <address@hidden>
>  > An: Rolf Fabian <address@hidden>
>  > CC: address@hidden
>  > Betreff: Re: left-right (up-down) parsing may lead to over/underflow of 
> prod
>
>
>
>  > Rolf Fabian 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.
>  > >
>  > > Rolf
>  > >
>  > >
>  > > -----
>  > > Rolf Fabian
>  > > <r dot fabian at jacobs-university dot de>
>  > >
>  > >
>  > This is consistent with MatlabR2007b. See
>  >
>  > >> x = [1e150, 1e160, 1e150, 1e160, 1e-310, 1e-310]
>  >
>  > x =
>  >
>  >   1.0e+160 *
>  >
>  >     0.0000    1.0000    0.0000    1.0000    0.0000    0.0000
>  >
>  > >> prod(x)
>  >
>  > ans =
>  >
>  >    Inf
>  >
>  > D.
>  >
>  > --
>  > David Bateman                                address@hidden
>  > Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph)
>  > Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob)
>  > 91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax)
>  >
>  > The information contained in this communication has been classified as:
>  >
>  > [x] General Business Information
>  > [ ] Motorola Internal Use Only
>  > [ ] Motorola Confidential Proprietary
>
>
>  What is actually your intention concerning the development of Octave ?
>
>  It seems to me that perfect replication of MatLab's bugs which even severly 
> contradict basic laws of algebra have higher precedence for
>  you than getting correct results.
>

It's not quite right to call this a bug, IMO. There is an obvious
speed-robustness tradeoff here, Matlab and Octave are just favoring
speed. Many such tradeoffs exist.
I basically share David's and Shai's view that it's OK for the builtin
`prod' to favor speed; still, it's desirable to address this in some
way.
After a little playing around I've come up with the attached
implementation of `safeprod' which seems to address the problem. Can
you check it out?
After possible improvements & fixes, the options are
1. commit the function to OctaveForge's general package
2. make it part of Octave
3. rewrite the algorithm in C++ and replace the builtin `prod'. Note
that this won't slow down too many cases, since the log2-based
algorithm (can be done using C's frexp) will only be used if the
"normal" reduction ends up with 0, Inf or NaN.

What do you think?

>  Rolf
>
>
>
>  --
>  GMX startet ShortView.de. Hier findest Du Leute mit Deinen Interessen!
>  Jetzt dabei sein: http://www.shortview.de/address@hidden
>



-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

safeprod.m
Description: Binary data