[jira] [Commented] (MATH-1558) MidpointIntegrator doesn't implement the midpoint method correctly

["Sam Ritchie (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/MATH-1558?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17215560#comment-17215560 ] 

Sam Ritchie commented on MATH-1558:
-----------------------------------

Ugh, sorry, the code formatting looked great in the editor but now looks goofy here on the issue.

I should note that the reason it was probably implemented this way is that doubling `n` _does_ work for the trapezoid method, if you want to reuse points:

{{n = 2: x-------x-------x}}
{{n = 4: x---x---x---x—--x}}

 

And TrapezoidIntegrator does implement this correctly. It looks like the MidpointIntegrator code was adapted from the TrapezoidIntegrator code without realizing the difference.

> 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
>
> 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:
> {{n = 1 |--------x--------|}}
>  {{n = 3 |--x--|--x--|--x–-|}}
> 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:
> {{n = 1 |-------x-------|}}
>  {{n = 2 |---x---|---x---|}}
>  {{n = 4 |-x-|-x-|-x-|-x-|}}
> The implementation of "stage" in MidpointIntegrator is combining them like this:
> {{n = 1 |-------x-------|}}
>  {{n = 2 |---x---x---x---|}}
>  {{n = 4 |-x-x-x-x-x-x-x-|}}
> 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:
> {{n = 2 x-------x--------}}
>  {{n = 4 x---x---x---x----}}
>  {{n = 8 x-x-x-x-x-x-x-x--}}
> Here's the code from "stage" implementing this:
> {{{{// 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);}}}}
> 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)