Numerical integration in Python with adaptive quadrature of vectorized function

Question

I'm looking for a super duper numerical quadrature function. It should have the following three properties:

• Adaptive - it automatically adjusts the density of sampling points to fit the integrand. This is absolutely necessary because my integrand is very nonuniform and expensive to compute.
• Vectorized - it calls the integrand on lists of sample points rather than one point at a time, for efficiency.
• Able to handle vector-valued functions - all components of the vector-valued integrand are computed at the same time for no additional cost, so it makes no sense to integrate all the components separately.

In addition, it should be:

• 2D - the integral I want to compute is a double integral over a planar region, and I want to be able to specify an overall (relative) tolerance for the whole integral and have it manage the error budget appropriately.

Does anybody know of a library that has such a function? Even two or three of the four properties would be better than nothing.

I'm using Python and SciPy, so if it already works with Python that's a bonus. (But I'm also able to write glue code to let it call my integrand if necessary.)

Answer 1

I just implemented vectorized adaptive quadrature for 1D and 2D domains in quadpy. All you need to provide is a triangulation of your domain and the function you want to integrate. It may be vector-valued.

Install quadpy with

pip3 install quadpy

and run

import numpy
import quadpy


triangles = numpy.array(
    [[[0.0, 0.0], [1.0, 0.0]], [[1.0, 0.0], [1.0, 1.0]], [[0.0, 1.0], [0.0, 1.0]]]
)

val, error_estimate = quadpy.triangle.integrate_adaptive(
    lambda x: [numpy.sin(x[0]), numpy.exp(x[0])], triangles, 1.0e-10
)

print(val)
print(error_estimate)

This gives

[ 0.45969769  1.71828183]
[  7.10494337e-12   3.68776277e-11]

Answer 2

The quadrature function in scipy.integrate satisfies the first two requirements of what you are looking for. The similar romberg function uses a different method.

Other functions only satisfy one of the requirements:

• The similarly-named quad function does adaptive quadrature, but only supports a function with a scalar argument. You can pass it a ctypes function for increased performance, but normal Python functions will be very slow.
• The simps function and related sampling methods can be passed a vector of (typically evenly-spaced) samples, but aren't adaptive.

The third requirement you listed (simultaneous integral of a vector-valued function) is a bit esoteric and conflicts with the ability to accept a vectorized function in the first place (the function argument would have to take a matrix!) Similarly, the ability to compute a double integral would complicate the specification of the function significantly.

In most cases, the quadrature function would be the way to go.

Answer 3

I used this library, it does everything you want, except it is written in C. But it also has an R interface, so maybe you can call R from Python (that is possible).

http://ab-initio.mit.edu/wiki/index.php/Cubature_(Multi-dimensional_integration)

Or, you can call the library using ctypes (it is not straight forward, but it is doable).