Adding an alternative integrator to deal with fast oscillations in long memory-length numerical integrations for CustomSD

[emile-cochin]

Motivation

I'm opening this issue following a discussion about potential problems encountered when computing process tensors with large memory lengths. The problem would arise from the numerical integration in CustomSD.eta_function due to the default integrator of scipy.integrate.quad not being able to deal with the fast oscillations of the integrand. Most of the times, oqupy manages to issue a warning to inform that there are potential errors in the computations but in niche cases, the errors are uncaught. I think this is because the influence functional computations using correlation_2d_integral involve (discrete) second order differences of eta_function which are very sensitive to tiny errors in eta_function.

The reason why I open this issue is to possibly contribute by adding an alternative to the current integrator to choose from if the current method doesn't yield satisfying results. It would rely on the weighted sin/cos integrators of scipy.integrate.quad.
Also, in order to diagnose problems with numerical integration, I believe that there is the helpers function plot_correlations_with_parameters for now. Unfortunately, sometimes the errors in the results of correlation_2d_integral don't quite appear from the bare correlations because of the sensitivity of the first/second order differences. It could maybe be useful to also add a function to plot the results for the 2d correlations as well.

Idea for the implementation

Just a quick remark first, I'm not sure that the documentation for eta_function is consistent with what it's computing. Shouldn't it be something like:

$$\eta(\tau) = \displaystyle\int_0^\infty \frac{J(\omega)}{\omega^2} \left[ (\cos(\omega\tau)-1)\coth\left(\frac{\omega}{2T}\right) - i(\sin(\omega\tau) - \omega\tau) \right]\mathrm{d}\omega$$

The current way this is computed is by splitting the integration interval into $[0, \omega_c]$ and $[\omega_c, +\infty]$ (for soft cutoffs). The method I tried which seemed to give convincing results when $\tau$ is larger than a few $1/\omega_c$ (if $\tau$ is too small, the integrand doesn't oscillate and using a weighted integrator doesn't work well) is the following:

• Further cut the first interval into $[0, \tilde{\omega}]$ and $[\tilde{\omega}, \omega_c]$
• Perform usual integration on the first interval $[0, \tilde{\omega}]$ because isolating the $\cos$ and $\sin$ parts of the integral isn't possible due to the integrands being divergent on their own. The choice of $\tilde{\omega}$ should be made such that there aren't too many oscillations in this interval/or at least that the default integrator is still able to integrate this part without error.
• For the two other intervals $[\tilde{\omega}, \omega_c]$ and $[\omega_c, +\infty]$, is it possible to split the integral into 4 parts:

$$\eta(\tau) = \displaystyle\int_a^b \frac{J(\omega)}{\omega^2}\coth\left(\frac{\omega}{2T}\right) \cos(\omega\tau)\mathrm{d}\omega - i \displaystyle\int_a^b \frac{J(\omega)}{\omega^2} \sin(\omega\tau)\mathrm{d}\omega - \displaystyle\int_a^b \frac{J(\omega)}{\omega^2} \coth\left(\frac{\omega}{2T}\right)\mathrm{d}\omega + i\tau\displaystyle\int_a^b \frac{J(\omega)}{\omega}\mathrm{d}\omega$$

The first two can be computed with the 'cos'/'sin' weighted integrals of quad. The other two aren't oscillating so they are computed with the default integrator. They also don't necessarily need to be recomputed if one doesn't change $\tilde{\omega}$ since they are independent of $\tau$. Right now, I'm not completely sure whether $\tilde{\omega}$ should be fixed by the user or estimated for each $\tau$ as something like $c/\tau$ where $c$ would be a fixed constant. For now, just fixing it at $\omega_c/2$ seemed to work for most of the spectral densities that I've tested but I don't know if it is desirable to have random constants sitting in the code like this.

Preliminary results

Here are a few examples for a some common spectral densities (which were chosen with an exponential cutoff which is where the usual integrator struggles most). The log-scale is a bit misleading but in most of these cases, there are still a few places where the errors are large enough to have an impact on the process tensor computations. Just as a side note, for the Ohmic case and the particular choice of parameters, the integration warning didn't trigger on the first large peak.

By a few days I'll commit the cleaned code. Let me know if you think this could be interesting to add, and if you have ideas of how to deal with choosing the ideal $\tilde{\omega}$.

Best,

[piperfw]

Hi @emile-cochin, thanks for the thorough work here - this seems like a very well posed issue.

Good spot on the eta_function docs, the current docstring is certainly incorrect (it seems to just be a copy of the autocorrelation function). Feel free to add that to your PR if you haven't already.

I'm away travelling but will look at the commited code soon. Your logic sounds good. Here are my main comments:

1. The choice of $\tilde{\omega}$. I agree we don't want a hidden magic number $c/\tau$. However, it may be reasonable to specify $\tilde{\omega}$ e.g. as an optional paramter or configuration variable that has a default value of $\omega_c/2$. The other possibility is to find a way to automate its estimation. Which one is the correct choice I think is a matter of whichever is more reliable. I.e., if $\omega_c/2$ really does work well for many cases, then probably using it as a default value but providing the user a (warning and) way to change it is easiest. Otherwise, a quick thought on automation is whether root finding can be used to calculate the number of oscillations reliably.

2. The results look really good. Regarding the side note: to be clear, are you talking about the scipy integration warning (as opposed to something one of us has written in OQuPy)? I wouldn't be surprised if that is unreliable. I'm not suggesting we try to develop an improved warning scheme over scipy, but maybe there are quad options e.g. epsrel, points that, if set appropiately, make the procedure more reliable? I'm not sure that would be easy to determine.

[emile-cochin]

Hi, thanks for the quick reply!

I've added a corrected docstring already for eta_function, let me know if it seems fine.

About the side note, I was indeed talking about the scipy integration warning. I'm not sure about whether there is an easy way to make the procedure more reliable but after trying many other parameters, this seemed to be the only unfortunate event where the warning wasn't triggered. Maybe there isn't much to do (other than keeping in mind that it's a good idea to verify the results from the numerical integration when trying long memory lengths).

For the choice of $\tilde{\omega}$, there is no real need for root finding since the oscillations are due to pure sine/cosine term $\cos(\omega\tau)$ or $\sin(\omega\tau)$. The number of oscillations on the interval $[0, \tilde{\omega}]$ is thus $n=\tilde{\omega}\tau/(2\pi)$. I think the clearest way of letting the user play with this $\tilde{\omega}$ parameter is to add a parameter num_oscillations to specify the number of oscillations kept for the first integral either as a constant or as a function of $\tau$. I've set $n=3$ as default value which generally gave very good results (although they are not very sensitive to the particular choice). The docs for these parameters is maybe a bit long but I didn't know how to make it short without referring to this particular GitHub issue.

All of this should now be committed.

[piperfw]

Hi emile, thanks for clarifying these points. Both sound fine to me. I'm a little time limited for the next two or three weeks but have made a note to review the code as soon as it can (or Gerald might before me). In the meantime, I added a couple of suggestions on the alt_integrator docstrings. I liked the None setting to 'inherit' the property (even though it seemed a bit weird to have the same parameter allowed different values between functions at first).

I don't think the length is the problem, and further referencing Issue #150 is not a bad idea - it may save 1 or 2 people a good deal of time trying to understand motivations for the parameter choice.

Other than that, have you thought about adding a short section/addition to the parameters tutorial on this?

[gefux]

Hi @emile-cochin,
Thank you very much for flagging & analysing these integration errors and delivering a solution along side it!
Concerning the codedesign I'll make a few suggestions now in the PR #151 thread.