commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Ajo Fod <ajo....@gmail.com>
Subject Re: [math] On MATH-995: Problems with LegendreGaussQuadrature class.
Date Tue, 02 Jul 2013 16:51:28 GMT
Konstantin,

In essence, we're both rooting for Adaptive Quadratures. IMHO
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 read. I
> 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.
>  http://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula
>
>
> 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]
> >
> > Please let me know if this doesn't answer your question about the
> > 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
> > > >> is necessary to be more informed about similar functionality 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 quadrature
> > > >>>> 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.]
> > > >>>>
> > > >>>> Additionally, your patch also mixes two concepts: Adaptive
> > > >>>> 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,
> please
> > > >>>>> 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".
> > > >>>>>>>>>
> > > >>>>>>>>>
> > > >>>>>>>>> Then please read all the messages 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