# commons-dev mailing list archives

##### Site index · List index
Message view
Top
From Konstantin Berlin <kber...@gmail.com>
Subject Re: [math] On MATH-995: Problems with LegendreGaussQuadrature class.
Date Tue, 02 Jul 2013 16:54:12 GMT
```IterateiveLegendreGaussIntegrator should be replaced by adaptive
guass-kronrod. The simpson and trapezoidal methods should be replaced by

On Tue, Jul 2, 2013 at 12:51 PM, Ajo Fod <ajo.fod@gmail.com> wrote:

> Konstantin,
>
> AdaptiveQuadrature (or some form of it that passes all Apache standards)
> must replace IterateiveLegendreGaussIntegrator as soon as possible before
> somebody else gets hit by a report of convergence to the wrong answer.
> **
>
> Cheers,
> -Ajo
>
>
>
>
> On Tue, Jul 2, 2013 at 9:33 AM, Konstantin Berlin <kberlin@gmail.com>
> wrote:
>
> > At first it seems you are still compute redundant points. See my example
> > that I posted, where I propagate 3 functional values not two.
> >
> > In regards to improvement. I am not an expert of different integration
> > strategies but:
> > The concept of adaptive quadrature is separate from how you integrate the
> > subinterval. In your implementation (and mine) we use Simpson's rule.
> This
> > is why putting comments on where your actual formulas come from is
> > important for public code.
> >         double Q1 = delta / 6 * (va + 4 * vm + vb); //this is simpson's
> > rule of integration
> > It does not make adaptive integration better or worse than Legendre-Gauss
> > quadrature, rather they solve different problems.
> > 2) The error tolerance should be an input to the quadrature, since
> > tolerance is a problem dependent value. This values should be subdivided
> > for each interval, such that the total sum would equal to the overall
> > desired tolerance. Look at the code I provided. In any case, like Gilles
> > said, the setup should be in the same spirit as other quadrature methods
> in
> > commons.
> > 3) Looking at Commons numerical integrations it is a little hard to
> > would have though that the basic quadratures, like trapezoidal or
> simpson's
> > rule, are adaptive already. I think its rare to use these basic formulas
> > without an adaptive quadrature. I would also think that Guass Kronrod
> would
> > be a better way of going about the integration than what is there now for
> > LGQ.
> >
> >
> > On Mon, Jul 1, 2013 at 11:37 PM, Ajo Fod <ajo.fod@gmail.com> wrote:
> >
> > > Hi Konstantin,
> > >
> > > Thanks for pointing out the inefficiency in AQ. I just improved the
> > > efficiency of AQ to 1.41x that of LGQ (up from 1.05x) - measured in
> > digits
> > > of accuracy per evaluation for integral of normal with sigma 1000 in
> > range
> > > [-5000, 5000]
> > >
> > > discussion:
> > >
> > > In essence Gilles thinks there is no problem with LGQ because it
> > integrates
> > > the low frequency functions in his unit tests accurately. He thinks
> that
> > > the problem is with the function I provided because it has "numerical
> > > instabilities" while I think it is the high frequency nature of the
> > > function that LGQ can't handle because it divides intervals
> > > indiscriminately. So, it is not clear to me how Gilles explains why AQ
> > > converges to the right answer in the presence of these "numerical
> > > instabilities" ... after all, LGQ and AQ are being passed the same
> > function
> > > and identical limits.
> > >
> > > This problem should appear with any function that has sufficiently high
> > > frequency components. It would be better if LGQ threw an exception when
> > it
> > > encounters a high frequency function. Instead, it "converges"
> confidently
> > > to the wrong answer. I personally think that AQ will fail at some point
> > at
> > > a high enough frequency ... but that will be well beyond the point at
> > which
> > > LGQ fails.
> > >
> > > I've avoided the complex method of fetching weights used in the current
> > > schemes because the  improvement in efficiency arises from the adaptive
> > > nature of the AQ method. You may notice that I'm using the weights that
> > you
> > > use in your code. I think Gilles requires that I use the weight
> > generation
> > > scheme he has worked with in the codebase in order to consider the code
> > > usable in Apache MATH.
> > >
> > > In summary, I feel the accuracy and versatility of AQ are being ignored
> > in
> > > favor of the familiarity of LGQ in the apache codebase. If there are
> > tests
> > > that AQ fails, I'll update my opinion.
> > >
> > > Cheers,
> > > Ajo
> > >
> > >
> > > On Mon, Jul 1, 2013 at 6:49 PM, Konstantin Berlin <kberlin@gmail.com>
> > > wrote:
> > >
> > > > I am not understanding the discussion here. Adaptive integration is
> > > > designed for functions that have different behavior in different
> > > > regions. Some regions are smoother, some have higher frequeniesy. How
> > > > you integrate a divided region, Simpson rule or whatever is a
> separate
> > > > question.
> > > >
> > > > Adaptive integration saves of function evaluations by avoiding large
> > > > series approximation in smooth regions. Nothing to do with how you
> > > > compute the subdivided regions.
> > > >
> > > > On Jul 1, 2013, at 8:45 PM,  Fod <ajo.fod@gmail.com> wrote:
> > > >
> > > > > Hi Gilles,
> > > > >
> > > > > Your accuracy concern made me wonder. So, I dropped the
> > > > > AdaptiveQuadrature.EPS to 1e-2 from 1e-9 in the code and ran the
> test
> > > in
> > > > > the patch.
> > > > >
> > > > > I computed the log of the error per evaluation ...i.e a measure of
> > the
> > > > > efficiency of the algorithm.
> > > > > And wait for it ... AQ beats LGQ by about 5% for the particular
> > > > formulation
> > > > > of the problem.
> > > > >
> > > > > Your request to use the Math classes falls under "coding style"
> IMHO.
> > > If
> > > > it
> > > > > doesn't satisfy your standards, feel free to modify. I'm happy with
> > it.
> > > > > Although as far as accuracy and convergence goes, I'd use AQ
> always.
> > > > >
> > > > >
> > > > > Got to compare apples to apples Gilles !
> > > > >
> > > > > Cheers,
> > > > > -Ajo
> > > > >
> > > > >
> > > > >
> > > > >
> > > > > On Mon, Jul 1, 2013 at 4:16 PM, Gilles <
> gilles@harfang.homelinux.org
> > >
> > > > wrote:
> > > > >
> > > > >> Hi.
> > > > >>
> > > > >>
> > > > >> On Mon, 1 Jul 2013 10:50:19 -0700, Ajo Fod wrote:
> > > > >>
> > > > >>> If you wanted to use the Math 3 codebase in AdaptiveQuadrature,
> > you'd
> > > > >>> compute the calculations of Q1 and Q2 with something else.
I'm
> not
> > > > >>> entirely
> > > > >>> familiar with the apache Math codebase [...]
> > > > >>
> > > > >> You could file a "wish" request as a Commons Math's user. Then,
if
> > > > >> and when some regular contributor finds some time, he will try
to
> > > > >> implement the functionality.
> > > > >>
> > > > >> However, when you provide a patch for inclusion in the codebase,
> it
that
> > > > >> would already exist in Commons Math, so that the contribution
can
> be
> > > > >> merged gracefully (i.e. with "minor" changes which committers
will
> > > > >> happily perform for you).
> > > > >> You are welcome to ask questions in order to be able to
> contribute.
> > > > >> As I tried to explain in more than one way, modifying your code
is
> > > > >> far from being trivial. If the committer has to figure out how
to
> > > > >> change/adapt/comment a significant part of the contribution,
it
> > > > >> ends up being easier to implement the feature from scratch!
> > > > >>
> > > > >>
> > > > >>
> > > > >>> Each of the tests in the patch is integrating a
> UnivariateFunction
> > in
> > > > >>> [-1,1]. Infinity.wrap(fn) just provides that UnivariateFunction.
> > [In
> > > > the
> > > > >>> patches for MATH-995 the InfiniteIntegral was replaced by
> > > > Infinity.wrap()
> > > > >>> ]. So, if you are saying that the intent of
> > > > >>> IterativeLegendreGaussIntegrat**or (refered to as LGQ) was
not to
> > > > >>> integrate
> > > > >>> this kind of UnivariateFunction in [-1,1], ... what kind
of
> > > univariate
> > > > >>> function would that be?
> > > > >>
> > > > >> Again, it is not just any UnivariateFunction, it is a function
> that
> > > > >> maps the [-inf, +inf] interval into [-1, 1].
> > > > >> It seems that the Gauss-Legendre quadrature is not appropriate
for
> > > > >> this. This is probably because the sample integration points
do
> not
> > > > >> cover the _whole_ interval: for the 10-point rule, the first
point
> > > > >> is at
> > > > >> -0.9739065285171717
> > > > >> and the last point is at
> > > > >>  0.9994069665572084
> > > > >> The interval in the original variable is thus [-18.908, 842.872].
> > This
> > > > >> is far from adequate for integrating a Gaussian function with
> > > > sigma=1000.
> > > > >> [And, as Phil pointed out from the outset, I suspect that the
> change
> > > of
> > > > >> variable also introduces numerical errors since the result becomes
> > > worse
> > > > >> when increasing the number of sample points. Increasing the
> > requested
> > > > >> precision leads to a prohibitive increase of the number of
> > > evaluations,
> > > > >> without improvement of the accuracy. In itself it is not
> sufficient
> > to
> > > > >> indicate a bug of "IterativeGaussLegendre"; it could simply be
a
> > > > >> limitation inherent to the algorithm.]
> > > > >>
> > > > >>
> > > > >> If it is indeed supposed to do the integration,
> > > > >>> then AQ clearly does a better job.
> > > > >>
> > > > >> Adaptive methods are certainly useful, but we need examples where
> > > > >> its usage is appropriate. It is _not_ indicated for the improper
> > > > >> integral of a Gaussian (even though it indeed performs better
than
> > > > >> Gauss-Legendre).
> > > > >>
> > > > >>
> > > > >>
> > > > >>> So, why does LGQ fail here? It is probably that the Adaptive
> > division
> > > > of
> > > > >>> the integration domain (as opposed to the uniform division
with
> > LGQ)
> > > > gives
> > > > >>> AQ the critical edge. The test you have for LGQ so far are
pretty
> > > well
> > > > >>> behaved.
> > > > >>
> > > > >> Cf. above.
> > > > >>
> > > > >> AFAIU, the problem reported by MATH-995 is not a bug in
> > > > >> "IterativeLegendreIntegrator": it correctly integrates a Gaussian
> > > > >> with a large sigma _if_ the integration interval is "large enough"
> > > > >> (cf. unit test referred to in my comment to MATH-995).
> > > > >>
> > > > >> Unless someone can point to something I'm missing in MATH-995,
> I'll
> > > > >> close that issue.
> > > > >>
> > > > >>
> > > > >>
> > > > >>> Summary: I'm demonstrating a clear bug/inefficiency with
LGQ
> > > > >>
> > > > >> I don't agree with that statement.
> > > > >> "IterativeGaussLegendre" produces the correct answer (at 1e-6
> > > > >> accuracy) in less than 60 function evaluations. To achieve the
> same,
> > > > >> your code needs 995 evaluations.
> > > > >>
> > > > >>
> > > > >> and providing
> > > > >>> you with an alternative that is more accurate.
> > > > >>
> > > > >> Cf. above (and my previous post), about how to contribute to
> > > > >> Commons Math.
> > > > >> Please open a new feature request.
> > > > >>
> > > > >>
> > > > >> Regards,
> > > > >> Gilles
> > > > >>
> > > > >>
> > > > >>> On Mon, Jul 1, 2013 at 8:22 AM, Gilles <
> > gilles@harfang.homelinux.org
> > > >
> > > > >>> wrote:
> > > > >>>
> > > > >>> Hi.
> > > > >>>>
> > > > >>>>
> > > > >>>>
> > > > >>>> I just noticed your request to write the algorithm along
the
> lines
> > > of
> > > > >>>>> the
> > > > >>>>> wikipedia article.
> > > > >>>>>
> > > > >>>>> The only major difference between my code and the
article on
> > > > Wikipedia
> > > > >>>>> is
> > > > >>>>> that I found it necessary to move the recursive stack
in into a
> > > data
> > > > >>>>> structure to avoid a StackOverflowException when
the non
> > polynomial
> > > > >>>>> curvature is concentrated in a corner of the domain
of
> > integration.
> > > > >>>>> Notice
> > > > >>>>> that the Stack objects stores a Stack of limits of
integration.
> > > > >>>> There is a misunderstanding: I'm referring to the "high-level"
> > > > >>>> description of the algorithm that is the separation of
concerns
> > > > >>>> between the quadrature method and the adaptive process.
Your
> code
> > > > >>>> mixes the two. Moreover, it does not reuse any of the
> > > > >>>> schemes already implemented in CM, but implements a (new?)
one
> > > > >>>> without any reference or comments.
> > > > >>>> [And this is even without delving into remarks concerning
the
> > > > >>>> code structure itself.]
> > > > >>>>
> > > > >>>> quadrature vs improper integral (which is also MATH-994);
it is
> > > > >>>> hard to follow what problem this issue is supposed to
point to,
> > > > >>>> and how the patch solves it. Indeed your unit tests shows
a
> > > > >>>> problem with improper integrals which the class
> > > > >>>> "****IterativeGaussLegendreIntegrat****or" is _not_ meant
to
> > > > handle.[1]
> > > > >>>>
> > > > >>>>
> > > > >>>> To be clear, hopefully, you are demonstrating a problem
that
> > > > >>>> occurs when combining Commons Math code with code which
you
> > > > >>>> created.
> > > > >>>> The first step is to create a unit test demonstrating
whether
> > > > >>>> an issue exists with "****IterativeGaussLegendreIntegrat****or"
> > code
> > > > >>>>
> > > > >>>> only (i.e. without relying on your "InfiniteIntegral"
class).[1]
> > > > >>>> If no independent issue exist, then MATH-995 should be
replaced
> > > > >>>> by an appropriate feature request.
> > > > >>>> Also, it would certainly be helpful to pinpoint the reason
why
> > > > >>>> the combination of "****IterativeGaussLegendreIntegrat****or"
> and
> > > > >>>>
> > > > >>>> "InfiniteIntegral" is not legitimate (if that's the case).
> > > > >>>>
> > > > >>>>
> > > > >>>> Regards,
> > > > >>>> Gilles
> > > > >>>>
> > > > >>>> [1] Cf. also my latest comment on the MATH-995 page.
> > > > >>>>
> > > > >>>>
> > > > >>>>
> > > > >>>> Cheers,
> > > > >>>>> Ajo.
> > > > >>>>>
> > > > >>>>>
> > > > >>>>> On Fri, Jun 28, 2013 at 11:07 AM, Ajo Fod <ajo.fod@gmail.com>
> > > wrote:
> > > > >>>>>
> > > > >>>>> BTW, it is possible that I'm not using LGQ correctly.
If so,
> > > > >>>>> show
> > > > >>>>>
> > > > >>>>>> how to pass the tests I've added. I'd much rather
use
> something
> > > > that is
> > > > >>>>>> better tested than my personal code.
> > > > >>>>>>
> > > > >>>>>> -Ajo.
> > > > >>>>>>
> > > > >>>>>>
> > > > >>>>>> On Fri, Jun 28, 2013 at 11:04 AM, Ajo Fod <ajo.fod@gmail.com>
> > > > wrote:
> > > > >>>>>>
> > > > >>>>>> I just posted a patch on this issue. Feel free
to edit as
> > > necessary
> > > > to
> > > > >>>>>>
> > > > >>>>>>> match your standards. There is a clear issue
with LGQ.
> > > > >>>>>>>
> > > > >>>>>>> Cheers,
> > > > >>>>>>> Ajo.
> > > > >>>>>>>
> > > > >>>>>>>
> > > > >>>>>>> On Fri, Jun 28, 2013 at 10:54 AM, Gilles
<
> > > > >>>>>>> gilles@harfang.homelinux.org>
> > > > >>>>>>> **wrote:
> > > > >>>>>>>
> > > > >>>>>>>
> > > > >>>>>>> Ted,
> > > > >>>>>>>
> > > > >>>>>>>>
> > > > >>>>>>>>
> > > > >>>>>>>>
> > > > >>>>>>>>  Did you read my other (rather more lengthy)
post?  Is that
> > > > >>>>>>>> "jumping"?
> > > > >>>>>>>>
> > > > >>>>>>>>>
> > > > >>>>>>>>>
> > > > >>>>>>>>>>
> > > > >>>>>>>>>> Yes.  You jumped on him rather
than helped him be
> > productive.
> > > > The
> > > > >>>>>>>>> general
> > > > >>>>>>>>> message is "we have something in
the works, don't bother us
> > > with
> > > > >>>>>>>>> your
> > > > >>>>>>>>> ideas".
> > > > >>>>>>>>>
> > > > >>>>>>>>>
pertaining to those
> issues
> > > more
> > > > >>>>>>>> carefully:
> > > > >>>>>>>> I never wrote such a thing (neither now
nor in the past).
> > > > >>>>>>>> I pointed to a potential problem in the
usage of the CM
> code.
> > > > >>>>>>>> I pointed (several times and in details)
to problems in
> > > candidate
> > > > >>>>>>>> contributions,
> > > > >>>>>>>> with arguments that go well beyond "bad
formatting".
> > > > >>>>>>>> I pointed out how we could improve the
functionality
> > _together_
> > > > (i.e.
> > > > >>>>>>>> by
> > > > >>>>>>>> using
> > > > >>>>>>>> what we have, instead of throwing it
out without even trying
> > to
> > > > >>>>>>>> figure
> > > > >>>>>>>> out how
> > > > >>>>>>>> good or bad it is).
> > > > >>>>>>>>
> > > > >>>>>>>> IMHO, these were all valid suggestions
to be productive in
> > > > helping CM
> > > > >>>>>>>> to
> > > > >>>>>>>> become
> > > > >>>>>>>> better, instead of merely larger. The
former indeed requires
> > > more
> > > > >>>>>>>> effort
> > > > >>>>>>>> than
> > > > >>>>>>>> the latter.
> > > > >>>>>>>>
> > > > >>>>>>>>
> > > > >>>>>>>>
> > > > >>>>>>>> Gilles
> > > > >>
> > > > >>
> > > > >>
> > > >
> > ------------------------------**------------------------------**---------
> > > > >> To unsubscribe, e-mail: dev-unsubscribe@commons.**apache.org<
> > > > dev-unsubscribe@commons.apache.org>
> > > > >> For additional commands, e-mail: dev-help@commons.apache.org
> > > > >>
> > > > >>
> > > >
> > > > ---------------------------------------------------------------------
> > > > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > > > For additional commands, e-mail: dev-help@commons.apache.org
> > > >
> > > >
> > >
> >
>

```
Mime
• Unnamed multipart/alternative (inline, None, 0 bytes)
View raw message