Only the first two sections of this long question are essential. The others are just for illustration.

Background

Advanced quadratures such as higher-degree composite Newton–Cotes, Gauß–Legendre, and Romberg seem to be mainly intended for cases where one can finely sample the function but not integrate analytically. However, for functions with structures finer than the sampling interval (see Appendix A for an example) or measurement noise, they cannot compete with simple approaches such as the midpoint or trapezoid rule (see Appendix B for a demonstration).

This is somewhat intuitive as, e.g., the composite Simpson rule essentially “discards” one quarter of the information by assigning it a lower weight. The only reason such quadratures are better for sufficiently boring functions is that properly handling border effects outweighs the effect of discarded information. From another point of view, it is intuitively clear to me that for functions with a fine structure or noise, samples that are remote from the borders of the integration domain must be almost equidistant and have almost the same weight (for a high number of samples). On the other hand, quadrature of such functions may benefit from a better handling of border effects (than for the midpoint method).

Question

Assume that I wish to numerically integrate noisy or fine-structured one-dimensional data.

The number of sampling points is fixed (due to function evaluation being costly), but I can freely place them. However, I (or the method) cannot place sampling points interactively, i.e., based on results from other sampling points. I also do not know potential problem regions beforehand. So, something like Gauß–Legendre (non-equidistant sampling points) is okay; adaptive quadrature is not since it requires interactively placed sampling points.

  • Have any methods going beyond the midpoint method been suggested for such a case?

  • Or: Is there any proof that the midpoint method is best under such conditions?

  • More generally: Is there any existing work on this problem?

Appendix A: Specific example of a fine-structured function

I wish to estimate $\int_0^1f(t)\, \mathrm{d}t$ for: $$ f(t) = \sum_{i=1}^{k} \frac{\sin(ω_i t-φ_i)}{ω_i},$$ with $φ_i∈ [0,2π]$ and $\log{ω_i} ∈ [1,1000]$. A typical function looks like this:

superimposed sines

I chose this function for the following properties:

  • It can be integrated analytically for a control result.
  • It has fine structure on a level that makes it impossible to capture all of it with the number of samples I am using ($<10^2$).
  • It is not dominated by its fine structure.

Appendix B: Benchmark

For completeness, here is a benchmark in Python:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(I here use the median to reduce the influence of outliers due to functions that have only high-frequency content. For the mean, the results are similar.)

The medians of the relative integration errors are:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

Note: After two months and one bounty without result, I posted this on MathOverflow.

  • Is this the kind of problem you really are interested in? In 1D, you can probably get good results rather quickly with most any method. – David Ketcheson Apr 27 at 6:18
  • "I have a fixed number of sampling points and can freely place them. However, I cannot place sampling points interactively, i.e., based on results from other sampling points." This restriction isn't clear to me. Am I allowed to put the nodes where an adaptive algorithm would put them, as long as I'm just really smart (instead of actually using the adaptive algorithm)? If I'm not allowed to be "really smart" about it, then what kind of node placements are actually allowed? – David Ketcheson Apr 27 at 6:31
  • @DavidKetcheson: Is this the kind of problem you really are interested in? – Yes, I am really interested in 1D. — In 1D, you can probably get good results rather quickly with most any method. – Remember that function evaluation may be costly. — then what kind of node placements are actually allowed? – I edited my question hoping to make it more clear. – Wrzlprmft Apr 27 at 8:02
  • Thanks that helps. To me, the question still seems vague. I think there is a simple and more precise question that would be more answerable. It would require defining a set of functions (that might depend on the allowed number of quadrature nodes) and a metric. Then you could ask if the midpoint method is optimal in that metric over that set of functions (where presumably the same set of nodes must be used for quadrature all functions). – David Ketcheson Apr 27 at 8:57
  • 1
    @DavidKetcheson: It would require defining a set of functions (that might depend on the allowed number of quadrature nodes) and a metric. – Given that I failed to find anything useful on this subject so far, I see no reason to impose such restrictions. Rather, with such restrictions, I would risk that I exclude some existing work (or easy proof) for slightly different conditions or assumptions. If there are any ways to capture the depicted scenario in definitions and similar for which a reference work or an easy proof exists, I am happy about that. – Wrzlprmft Apr 27 at 9:25

First of all, I think you misunderstand the concept of adaptive quadrature. Adaptive quadrature does not imply "interactively placing sample points". The whole idea behind adaptive quadrature is to devise a scheme that will integrate a certain function to a certain (estimated) absolute or relative error with as little function evaluations as possible.

A second remark: you write "The number of sampling points is fixed (due to function evaluation being costly), but I can freely place them". I think the idea should be that the number of sampling points (or function evaluations in quadrature terminology) should be as small as possible (i.e. not fixed).

So what is the idea behind adaptive quadrature as implemented in QUADPACK for example?

  1. The basic ingredient is a "nested" quadrature rule: this is a combination of two quadrature rules where one has a higher order (or accuracy) as the other. Why? Based on the difference between these rules, the algorithm can estimate the quadrature error (of course, the algorithm will use the most accurate one as the reference result). Examples could be the trapezoid rule with $2^{n}$ nodes and $2^{n+1}$ nodes. In the case of QUADPACK, the rules are Gauss-Kronrod rules. These are interpolatory quadrature rules that use a Gauss-Legendre quadrature rule of a certain order $N$ and an optimal extension of this rule. This means that a higher quadrature order can be obtained by re-using the Gauss-Legendre nodes (i.e. the costly function evaluations) with different weights and adding a number of extra nodes. In other words, the original Gauss-Legendre rule of order $N$ will integrate all polynomials of degree $2N-1$ exactly while the extended Gauss-Kronrod rule will integrate some higher order polynomial exactly. A classic rule is the G7K15 (7th order Gauss-Legendre with 15th order Gauss-Kronrod). The magic is that the 7 nodes of the Gauss-Legendre are a subset of the 15 nodes of the Gauss-Kronrod so with 15 function evaluations, I have a quadrature evaluation together with an error estimate!

  2. Next ingredient is a "divide-and-conquer" strategy. Suppose you let loose this G7K15 on your integrand and you observe a quadrature error that is according to your taste too large. QUADPACK will then subdivide the original interval in two equally spaced subintervals. And then it will re-evaluate the two subintegrals using the basic rule, G7K15. Now, the algorithm has a global error estimate (which should be hopefully lower than the first one) but also two local error estimates. It picks the interval with the largest error and divides this one in two. Two new integrals are estimated and the global error is updated. And so on until the global error is below your requested target or the maximum number of subdivisions has been surpassed.

So I challenge you to update your code above using the scipy.quad method. Maybe in case of an integrand with a lot of "fine structure" you might need to increase the maximum number of subdivisions (the limit option). You could also play with the epsabs and/or epsrel parameters.

However, if you only have experimental data, I see two possibilities.

  1. If you have the opportunity to select the measurement points, i.e. values of $t$, I would select them equidistanly and preferably as a power of $2$ so that you can apply a nested trapezoidal rule (and profit from Romberg extrapolation).
  2. If you have no means of choosing the nodes, i.e. the measurements come at random times, the best option in my opinion still is the trapezoid rule.
  • I think you misunderstand the concept of adaptive quadrature. – Your post completely agrees with my previous understanding of adaptive quadrature and it’s a clear match for how I defined interactively placing the sampling points (whether that’s an appropriate phrase or not). — you write […]. I think the idea should be that the number of sampling points […] should be as small as possible (i.e. not fixed). – If you have that luxury, sure, but experimental constraints may not be that benign. E.g., suppose you have to measure something simultaneously with a fixed number of expensive sensors. – Wrzlprmft Sep 4 at 23:08
  • My apologies. I misinterpreted "interactively" in your question. In my understanding "interactively" means intervention by the user not by an algorithm. I have added a paragraph in my answer on experimental data. Another approach would be to "filter" out the fine structure information, i.e. apply a Fourier transform and remove high-order frequencies with small amplitudes. Would that be an option? – GertVdE Sep 7 at 6:51
  • If you have the opportunity to select the measurement points […] – Equidistant points are what I need for midpoint, plain trapezoid, etc. anyway, so this is exactly what I did in my benchmark. Here, Romberg extrapolation does not yield any advantage. – Wrzlprmft Sep 7 at 9:58
  • Another approach would be to "filter" out the fine structure information […] Would that be an option? – In my example, I assume the fine structure is part of what I want to measure, I just do not happen to have sufficiently many samples to capture it completely. As for actual noise, there is no technical constraint that keeps me from filtering. However, the integral over the entire domain is already the ultimate low-pass filter, so I am skeptical that this can be improved without having noise with specific, benign, and known properties. – Wrzlprmft Sep 7 at 10:11
  • Is it truly stochastic? There must be some derived which are higher order stochastic integral approximations. – Chris Rackauckas Sep 7 at 17:23

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.