[MATH-1558] MidpointIntegrator doesn't implement the midpoint method correctly

• Type: Bug
• Status: Open
• Priority: Major
• Resolution: Unresolved
• Affects Version/s: 3.6.1
• Labels: integration midpoint
• Assignee: Unassigned
• Reporter: Sam Ritchie
• Created: 16/Oct/20 17:34
• Updated: 21/Oct/20 00:34

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.

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);

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!