[
https://issues.apache.org/jira/browse/MATH-1558?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17217982#comment-17217982
]
Sam Ritchie edited comment on MATH-1558 at 10/20/20, 10:36 PM:
---------------------------------------------------------------
Whoops, [~erans] somehow it escaped me that the project lived on GitHub. I've
just made a PR and added some tests verifying that we are saving roughly 1/3 of
evaluations vs a non-incremental approach. The rest of the tests pass (they
weren't making particularly strong claims about the results, just that
iterations didn't exceed the max).
I've also just sent an ICLA to [email protected], so we should be all set
there too.
was (Author: sritchie):
Whoops, [~erans] somehow it escaped me that the project lived on GitHub. I've
just made a PR and added some tests verifying that we are saving roughly 1/3 of
evaluations vs a non-incremental approach. The rest of the tests pass (they
weren't making particularly strong claims about the results, just that
iterations didn't exceed the max).
> MidpointIntegrator doesn't implement the midpoint method correctly
> ------------------------------------------------------------------
>
> Key: MATH-1558
> URL: https://issues.apache.org/jira/browse/MATH-1558
> Project: Commons Math
> Issue Type: Bug
> Affects Versions: 3.6.1
> Reporter: Sam Ritchie
> Priority: Major
> Labels: integration, midpoint
> Attachments: midpoint.patch
>
> Time Spent: 10m
> Remaining Estimate: 0h
>
> Hi all,
> I've been reading through the implementation of {{MidpointIntegrator}}, and
> I've discovered an issue with the algorithm as implemented.
> The midpoint method generates an estimate for the definite integral of {{f}}
> between {{a}} and {{b}} by
> * subdividing {{(a, b)}} into {{n}} intervals
> * estimating the area of each interval as a rectangle {{(b-a)/n}} wide and
> as tall as the midpoint of the interval - ie, {{((b-a)/n) * f((a + b) / 2)}}
> * adding up all estimates.
> The {{MidpointIntegrator}} implementation here sticks to n = powers of 2 for
> this reason, stated in the comments:
> ?? The interval is divided equally into 2^n sections rather than an arbitrary
> m sections because this configuration can best utilize the already computed
> values.??
>
> *Here is the* *problem:* the midpoint method can't reuse values if you split
> an interval in half. It can only reuse previous values if you split the
> interval into 3; ie, use 3^n sections, not 2^n.
> What's actually implemented in `MidpointIntegrator` is a left Riemann sum
> without the leftmost point. Or, identically, a right Riemann sum without the
> rightmost point: [https://en.wikipedia.org/wiki/Riemann_sum]
> This matters because the error of a (left, right) Riemann sum is proportional
> to 1/n, while the error of the midpoint method is proportional to 1/n^2... a
> big difference.
> h2. Explanation
> The idea behind the midpoint method's point reuse is that if you triple the
> number of integral slices, one of the midpoints with n slices will overlap
> with the n/3 slice evaluation:
> {noformat}
> n = 1 |--------x--------|
> n = 3 |--x--|--x--|--x–-|
> {noformat}
> You can incrementally update the n=1 estimate by
> * sampling the points at (1/6) and (5/6) of the way across the n=1 interval
> * adding them to the n=1 estimate
> * scaling this sum down by 3, to scale down the widths of the rectangles
> For values of n != 1, the same trick applies... just sample the 1/6, 5/6
> point for every slice in the n/3 evaluation.
> What happens when you try and scale from n => 2n? The old function
> evaluations all fall on the _boundaries_ between the new cells, and can't be
> reused:
> {noformat}
> n = 1 |-------x-------|
> n = 2 |---x---|---x---|
> n = 4 |-x-|-x-|-x-|-x-|
> {noformat}
> The implementation of "stage" in MidpointIntegrator is combining them like
> this:
> {noformat}
> n = 1 |-------x-------|
> n = 2 |---x---x---x---|
> n = 4 |-x-x-x-x-x-x-x-|
> {noformat}
> which is actually:
> * tripling the number of integration slices at each step, not doubling, and
> * moving the function evaluation points out of the midpoint.
> In fact, what "stage" is actually doing is calculating a left Riemann sum,
> just without the left-most point.
> Here are the evaluation points for a left Riemann sum for comparison:
> {noformat}
> n = 2 x-------x--------
> n = 4 x---x---x---x----
> n = 8 x-x-x-x-x-x-x-x--
> {noformat}
> Here's the code from "stage" implementing this:
> {code}
> // number of new points in this stage... 2^n points at this stage, so we
> // have 2^n-1 points.
> final long np = 1L << (n - 1);
> double sum = 0;
> // spacing between adjacent new points
> final double spacing = diffMaxMin / np;
> // the first new point}}}}
> double x = min + 0.5 * spacing;}}}}
> for (long i = 0; i < np; i++) {}}}}
> sum += computeObjectiveValue(x);
> x += spacing;}}}}
> }
> // add the new sum to previously calculated result
> return 0.5 * (previousStageResult + sum * spacing);
> {code}
> h2. Suggested Resolution
> To keep the exact same results, I think the only solution is to remove the
> incorrect incremental "stage"; or convert it to the algorithm that implements
> the correct incremental increase by 3, and then _don't_ call it by default.
> Everything does work wonderfully if you expand n by a factor of 3 each time.
> Press discusses this trick in Numerical Recipes, section 4.4 (p136):
> [http://phys.uri.edu/nigh/NumRec/bookfpdf/f4-4.pdf] and notes that the
> savings you get still make it more efficient than NOT reusing function
> evaluations and implementing a correct scale-by-2 each time.
> Happy to provide more information if I can be helpful!
>
--
This message was sent by Atlassian Jira
(v8.3.4#803005)
