[
https://issues.apache.org/jira/browse/MATH-1558?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel
]
Gilles Sadowski updated MATH-1558:
----------------------------------
Description:
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!
was:
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!
> 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:
> {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)
