I’m interesting in calculating an integral of the form
\int_{\mathbb{R}^n} d^n x f(x,\beta) \exp(-\beta H(x))
where H is a known an Hamiltonian, \beta is an inverse temperature and f is a function which I can only calculate numerically. For that reason, I don’t have much information about its behavior, but it seems at least smooth. For now, I’m working with n=4, but, in the future, I would be interested in going to higher dimensions.
For a certain f, I had success in truncating the integration region and using HCubature from Integrals.jl. The thing is that, as I increase \beta, the integrand becomes more concentrated and I’m forced to “update” the integration interval more or less by hand. The automatic change of variables that is used by Integrals.jl to handle infinite domains didn’t work for me, because it simply gives zero to my integral. For other f, I couldn’t make this truncation work.
I was also able to perform a Monte Carlo integration by sampling points from the density \exp(-\beta H(x)) . This works for all the f 's that I tested but I need a LOT more integration points (\sim 10^7) to get consistent results.
Would anyone have ideas on better ways to tackle this problem? Thanks!
Can you do a substitution to scale out the concentration of measure? Perhaps you can change your integration variable from x to u = x / \lambda where \lambda is the thermal de Broglie wavelength or some other temperature-dependent length scale in your Hamiltonian?
Problably you just need a rescaling by some units so that it samples at an appropriate scale. See the comments about lengthscale in the QuadGK example for infinite limits.
I’ve seen problems with multiple length scales (temperature included) where infinite limits aren’t possible due to limited datasets. In those cases, when truncating the domain can be done in an automated way it has worked well. Since you mentioned H is a known Hamiltonian, it might be possible to do this. For example, in 1d this means finding the interval in which the Boltzmann factor is larger than \epsilon>0. If for some function, G, we know G(x) \geq H(x), \forall x and G is easily invertible, then the domain on which \exp(-\beta H(x)) > \epsilon > 0 is a subset of the interior of the level set G^{-1}(\log(1/\epsilon)/\beta). This procedure works best if G is an asymptotically tight upper bound on H, or it could be H itself.