Thanks for the post and resources. I wasn’t sure if anyone would be interested in mentoring the following, but it’s a project I’m interested in working on related to my research. I’ll leave the details below, but will try and go to the Fortran meeting tomorrow.
Something I am interested in doing is creating a multidimensional integration routine, for the below form of integrals. If there were an interested mentor, perhaps this could be something to contribute as a Fortran project?
Problem
In general, the denominator may approach 0 when f does not; depending on the number of dimensions, this is in fact an integrable singularity. However, traditional integration methods take a large number of points to be able to resolve this region accurately. In addition, these functions f and g may only be known numerically on a given set of points (due to being difficult to compute), so adaptive mesh refinement may not be possible, or if it’s being used, it’s at a higher-level. This is the use case I’ve been thinking of.
In condensed matter physics, these are typically 2D or 3D integrals over a convex domain, and the closest thing they have to this is the so-called Tetrahedron Method, which is only able to compute the imaginary part of the above integrals, for z = i\epsilon; additionally, the tetrahedra meshes must be uniform. Additionally, they rely on f(x) = 1. Despite these restrictions, this so-called Tetrahedron method is amazing because of the massive reduction in sampling points needed.
Scuffed Example
As an example, here is an early comparison I made a few months ago of computing the above integrals in 3D with f=1, g(x) ~ cos(ax) + cos(by) + cos(cz), plotting I(w+i eps) for small epsilon.:
The naive method essentially is the rectangle method. We get basically equivalent accuracy between the two methods, although one mesh is 4^3 as small. [The asymmetry in the tetrahedron method result came from an accidental shift of the mesh that produced a systematic error] Additionally, in the limit epsilon → 0, the rectangle method becomes super noisy, whereas the tetrahedron does not.
Improvement/Project Idea
However, Kaprzyk in this set of papers:
3D algebraic (generalized) tetrahedron method
N-D algebraic (generalized) simplex method
outlines a procedure for an integration routine that works for non-uniform meshes of simplices, works for arbitrary f(x) and g(x) (including numerically generated ones), and generalized to N-dimensions. In the 3D case, he has hardcoded his idea for a given set of functions, and there is no code for the N-D case.
The idea is as such: you subdivide the domain into simplices, linearly interpolate f(x) and g(x) in each simplex, and use analytic integration formulae to evaluate the integrals. Each point in the domain gets a weight, and the most difficult part of the computation becomes evaluating these parts of the weights:
The above are the formulae for 3-D integrals, but there is a similar form (with bigger determinants and higher powers of z) for N-dimensions.
Pros/Advantages of the Method
- From the simple experiment above, getting the same accuracy for a 4^N smaller mesh seems great
- The BIGGEST strength is that the integrals w.r.t. z of I(z) above are also analytically obtained with a simple change to the weights computed above. What this means is you can also compute integrals of I(z) without needing to repeat the above for N_z mesh points.
- In other words, computing \int I(z) dz (and further integrals) is the same cost as computing I(z). This will be very useful for many-body condensed matter applications (the GW method, electron-phonon calculations, etc) and I’m sure would be generally applicable.
- The formalism is agnostic to the chosen mesh (uniform vs nonuniform, etc)
Cons of the Method
- The user has to have the numerator f(x) and denominator g(x) known separately.
- If the denominator doesn’t approach 0, then in general other integration schemes will be faster or more accurate.
Issues to solve for the project:
- Obtaining simplified formulae for when the quantities z_j → z_k, which may occur in some simplices, for high dimensions will require a fair amount of algebra. However, I believe I can simplify the process better than Kaprzyk did.
- Evaluating the above quickly. They in general simplify to (N-th recursive) finite differences [if all z_i are distinct) which ends up being fairly slow. A clever lookup table/memoization could massively speed up this part of the calculation.
- Interfacing with an N-dim Delaunay generator; Qhull appears to be the only library that does this for N-dimensions. Up to 8-D integrals might be the max that’s realistic however.
