Hello all! This issue is seeking thoughts on a core part of scipy.integrate. More than that, it is proposing an enhancement that I and others are willing and able to work towards.
Integrating a function $f(x)$ is a fundamental part of numerical computing, and doing so for a function $f(\mathbf{x})$ in $n$ dimensions is just as important as doing so in 1 dimension. SciPy currently has routines in its integrate submodule to do so, notably dblquad, tplquad, and nquad. These are wrappers that essentially do 1D iterated integrals using quad from the Fortran library QUADPACK. QUADPACK is very impressive software, but at the end of the day it is still written in Fortran and limited by the design of its time. The version in SciPy is from 1984.
Improving integrate generally is important, and I'd very much like feedback on a plan on how to do that for quadrature in $n$ dimensions. If what is written here is broadly acceptable to everyone, I and others would intend to carry out this work on multidimensional integration before the end of this summer. The inspiration here is really the adaptive quadrature in quad_vec, which I think went a long way to providing robust integration in SciPy that works well with arrays.
To be very specific about the issues, multidimensional integration in SciPy currently has these limitations:
It does not support array-valued integrands of arbitrary shape, only scalar functions. Therefore one cannot integrate a vector or matrix or tensor function whereas quad_vec can.
It does not support vectorization. If one wants to speed things up by sampling all the integration points at once, you can't.
It will never support other array APIs, as it is written in Fortran. So it can never go on a GPU.
It does not support custom integration rules. If one wants to use, say, multidimensional tanh-sinh quadrature for handling a singular point, you can't.
The proposal here is to provide a multidimensional integration function that is adaptive and similar to quad_vec. This is what is known as "cubature", we can call it cub or ndquad_vec or whatever. It will take almost the same arguments as quad_vec, except a, b, and points will all become $n$ dimensional. It will use the array API to enable support for GPU integration and beyond.
An additional feature that it will have, which quad_vecdoesn't have, is support for passing in a custom integration rule. This would allow the adaptive part to subdivide effectively, but it could then use the custom integration rule to calculate the integral in each subdivision. The built-in rules would be Gauss-Kronrod and Genz-Malik (popular for multidimensional integrands). But users could also pass in rules that effectively deal with singular (using tanh-sinh) or oscillatory integrands. QUADPACK actually has behaviour like this with its methods "qawoe", "qawse", and "qawce", but it is not fully custom. If successful, this behaviour could be added into quad_vec as well by simply making its existing quadrature argument also accept a callable.
This proposal would go a long way towards improving an existing capability within SciPy. It would also help a lot with removing Fortran from the codebase.
I know @mdhaber, among others, has been doing work on integrate. I'd very much appreciate his thoughts and, of course, everyone's.
vnmabus, jorenham, DavidMStraub and BearBearCodeslucascolley
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4