Integration of functions over geometries

[stevengj]

@mike.ingold

It looks like the implicit caveat to the “16-digit accuracy and in O(n) time” bumper sticker is that this is only precisely true for functions that can be exactly represented by a polynomial of order 2n-1 over the specified interval.

No, that’s not the case. For polynomials Gaussian quadrature is exact, but for analytic functions Gaussian quadrature converges exponentially fast with the number of quadrature points n, so you can indeed get to nearly machine precision (nearly 16 digits) rapidly, and if you need to go to large n then FastGaussQuadrature.jl computes the quadrature rule in O(n) time.

A downside of FastGaussQuadrature.jl is that it just gives you an n-point rule, it doesn’t tell you what n should be for a given tolerance and a given function, and the rule is not a “nested” rule — it doesn’t give you an error estimate. So you are only your own in finding the correct n for your problem, e.g. by repeatedly doubling n until it converges. (Repeatedly increasing the order is also called p-adaptive quadrature.)

QuadGK is h-adaptive, in that it repeatedly subdivides the integration interval until a tolerance is achieved, and it uses an embedded rule to get an error estimate. This is convenient, and it is also fairly robust if your integrand is badly behaved in some localized region because it will only refine the quadrature around that region.

If your function is analytic, then using high-order quadrature is theoretically more efficient, but on the other hand you rarely have to go to very high order in practice, and the time for the computation of the rule is less important if you are doing integration over and over for many different functions and/or domains. So for smooth functions you can simply pass order=n for a larger n to QuadGK.jl to use a higher-order rule (which is computed once and cached). The computation of the rule is O(n^2), so it will be much slower than FastGaussQuadrature.jl for large n, but since in practice you rarely need n larger than a few hundred it doesn’t matter too much in my experience.

The other big advantage of FastGaussQuadrature.jl is that it has built-in rules for weighted quadrature with many different weight functions, which allow you to efficiently handle known endpoint singularities. See the discussion of weighted Gauss rules in the QuadGK manual, and integrands with singularities and discontinuities.

PS. Although QuadGK is not asymptotically exponentially convergent for a fixed order, unlike repeatedly increasing the order, it often gives the illusion of exponential convergence up to a given precision. See also this post: (Insanely-) Higher Order Derivatives - #17 by stevengj

[mike.ingold]

It’s for exactly these reasons that I generally tend toward h-adaptive quadrature/cubature methods. These days I find myself writing a lot of tooling that will be used to run calculations/models, moreso than direct data-in-hand problems. That being the case, I find it hugely beneficial to have an interface where I can specify some acceptable atol or rtol and get a result within those specs reliably since I can’t know at time of writing what function/signal will even be integrated.

[mike.ingold]

Updates on LineIntegrals:

• I’ve re-implemented all of the existing integral methods using FGQ.
• All QuadGK.quadgk(f, <:Meshes.Geometry) methods remain defined. Again, this is type piracy so would probably need to be PR’d into another package.
• I’m planning to implement an internal function with in-place memory re-use for geometries with multiple segments.
• The README has a toy example integrating f(x,y) = abs(x+y) over a Bezier curve approximation of the unit circle using both techniques. The integrate methods default to n=100 for Gauss-Legendre nodes, and quadgk uses all default kwargs. Both are solved in similar time with similar allocations, but QuadGK.jl seems to be more accurate. Increasing n with integral causes its result to converge toward quadgk’s, but at the expense of a significant increase in allocations. Not a scientific test by any means; your mileage may vary.

[lxvm]

Hi,

I’m new to the thread, so I may have missed some details, and I wanted to share a package that allows integration over polytopes using nested integration as an algorithm, which is convenient for integration in low-dimensions (up to 3 or 4), and even ideal for nearly singular integrands. In two (and more) dimensions the nested integration algorithm evaluates the multivariate integral in terms of recursive 1d integrals, i.e:

\int_{a_x}^{b_x} dx\ \int_{a_y(x)}^{b_y(x)} dy\ f(x,y) = \int_{a_x}^{b_x} dx\ I_y(x), \quad I_y(x) = \int_{a_y(x)}^{b_y(x)} dy\ f(x,y)

What is nice about this algorithm is that it is compatible with any quadrature scheme and any geometry (polytopes, manifolds) as long as either a parametrization of the geometry is available, or optimization on that geometry is possible and fast compared to evaluation of the integrand.
The example below is with the package AutoBZCore.jl, which emulates the Integrals.jl interface, and it allows for integration over polytopes from Polyhedra.jl

julia> using AutoBZCore, Polyhedra

# make convex hull in 3d
julia> domain = load_limits(polyhedron(vrep(rand(10, 3))))
PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078]))

# setup integrand
julia> prob = AutoBZCore.IntegralProblem((x,p) -> sin(prod(x)), domain)
AutoBZCore.IntegralProblem{var"#5#6", PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}, AutoBZCore.NullParameters}(var"#5#6"(), PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078])), AutoBZCore.NullParameters())

# call quadgk recursively
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.019284903620528718, 1.7385351242394514e-18, true, -1)

julia> using FastGaussQuadrature

# use a product rule of 10-point Gauss-Legendre quadratures
julia> AutoBZCore.solve(prob, NestedQuad(QuadratureFunction(fun=gausslegendre, npt=10)))
AutoBZCore.IntegralSolution{Float64, Nothing}(0.01928490362052873, nothing, true, -1)

# pick quadgk for two variables and Gauss-Legendre for the last
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL(), QuadratureFunction(fun=gausslegendre, npt=10), QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.01928490362052872, 5.454892180317694e-19, true, -1)

My hope is to integrate this functionality with Integrals.jl (no pun intended) when I find the time, and the main idea of the implementation is in
IteratedIntegration.jl.

[mike.ingold]

Nice work! I’d actually been wondering earlier, among the discussions about the different integration packages, if the better interface wouldn’t be implement this functionality into Integrals.jl. The core functionality that I’m really looking for, personally, is the ability to integrate functions over some standardized 1D/2D/3D geometry definitions, ideally with an easy way to derive the geometry from a CAD mesh model.

[juliohm]

@mike.ingold

The core functionality that I’m really looking for, personally, is the ability to integrate functions over some standardized 1D/2D/3D geometry definitions, ideally with an easy way to derive the geometry from a CAD mesh model.

Same here @mike.ingold. The ability to integrate functions over CAD-derived geometries is really missing, and you are closing this gap

As soon as we finish the new coordinate system infrastructure in Meshes.jl, we will be able to address very practical problems such as computing areas and volumes of complicated domains, which are not made of polytopes.

[lxvm]

@juliohm Out of curiosity, does Meshes.jl implement elimination, i.e. finding a lower-dimensional hull from the intersection of a polytope and a hyperplane, and double description, i.e. computing half-spaces from vertices and vice-versa? Both of these are all you need for applying iterated integration to convex polytopes.

[juliohm]

@lxvm we don’t have those. It would be nice to review PRs with these additions

[aplavin]

Does “integration of functions over geometries” include 2d integration over arbitrarily defined regions, through automatic variable change?
Like, computing \int_a^b dx \int_{lo(x)}^{hi(x)} f(x, y) dy.

[juliohm]

@aplavin as discussed above, we are considering all sorts of integrations over geometries. People are working on different fronts to make these more readily available.

[lxvm]

Yes, that is the goal. In this direction, IteratedIntegration.jl defines an interface for that form of limits of integration and uses Polyhedra.jl to implement it for integration over polytopes. The package is still experimental in its algorithms, but that interface on integration domains will remain stable and could be extended to more geometries that are convex and support optimization for finding the extrema of the coordinate axes on that manifold.