[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: left-right (up-down) parsing may lead to over/underflow of prod
From: |
Jaroslav Hajek |
Subject: |
Re: left-right (up-down) parsing may lead to over/underflow of prod |
Date: |
Tue, 22 Apr 2008 21:53:46 +0200 |
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
Re: left-right (up-down) parsing may lead to over/underflow of prod, Shai Ayal, 2008/04/22