Recently, I tried to numerically integrate a function for multiple parameters. The integral is four dimensional, so I use CubatureJLh() (but I don’t mind using other methods - the documentation doesn’t elaborate a lot on the differences, so I just chose one that supports multidimensional integration and batching).
The parameters I’m passing to the integral affect both the integrand itself and the integration boundaries (of some of the dimensions), and I need to evaluate a large matrix of parameters (about 100x100 pairs of parameters).
As a first step, I ran the computation (for each value of the parameters) in a multithreaded manner. This was sufficient when integrating a simple function, but for another integrand that I needed, this is still very slow. Generally, my integral consists of a product of exponents of a (smooth) function.
Integrals.jl’s caching interface automates this process to reuse resources if an algorithm supports it and if the necessary types to build the cache can be inferred from prob.
My questions are:
Which methods, other than QuadJKL (which only supports one-dimensional integrals), supports caching?
More generally, when can I expect caching to help with the performance? So far, I’ve used numerical integration as a blackbox, so it’s not immediately clear to me how caching of the computation of the integral for one parameter set can help compute the same integral for other values of the parameters.
Other general tips for speeding up the computation for this kind of problem (multidimensional integral and a large matrix of parameters) would be appreciated.
I believe HCubatureJL, which is a pure-Julia rewrite of CubatureJLh, allows for some caching. (It mainly stores a heap for h-adaptive refinement.) It doesn’t support batching, but it does allow integration of arbitrary types, so you can efficiently batch with a vector-valued integrand using the StaticArrays.jl package and pass an arbitrary norm function for computing the batch error as you see fit (e.g. L_2 or L_inf).
Caching in integration packages is typically done in two ways: caching of quadrature rules, and caching of data structures used in domain decomposition. The former is used in p-adaptive quadrature where the degree of the quadrature rule is increased until convergence is reached, and the latter is used in h-adaptive quadrature where the same quadrature rule is applied to each piece in a tiling of the integration domain, subdividing the tiles until convergence is reached. p-adaptivity is ideal for smooth functions and h-adaptivity is ideal for functions with localized features.
Speeding up integrals is highly problem-dependent. You may find that black box algorithms work well enough for your problem, but by exploiting the structure of your integrand (e.g. periodicity in certain variables or known functional forms) you can likely build a specialized method that works better
The integral is four dimensional
In 4D both h-adaptive and p-adaptive quadrature will be fairly costly, but still more accurate than Monte Carlo methods. Specialized quadrature methods in 4d are rare.
Generally, my integral consists of a product of exponents of a (smooth) function.
This is an indication that p-adaptive quadrature may perform well for your problem. CubatureJLp uses Clenshaw-Curtis quadrature to do this, although you can easily form tensor product quadratures, e.g. using FastGaussQuadrature.jl, and write logic for p-adaptive refinement by comparing two quadratures with different degrees until convergence is reached.
There are specialized quadratures for functions with known singularity locations and types, functions that are highly oscillatory, and for integrals with a known weight function to name a few.