I need to numerically compute an integral of the following form: $$\int_0^1 \frac{1}{2\pi\sigma^2}\exp\left(-\frac{\|(q_0t^3 + q_1t^2 + q_2t + q_3) - a\|^2}{2\sigma^2}\right)\|3q_0t^2 + 2q_1t + q_2\|\,dt.$$ Here, $q_0,q_1,q_2,q_3,a$ are constant vectors in $\mathbb{R}^2$, and $\sigma$ is real. This comes from the line integral of a Gaussian over a cubic curve.

Obviously I can integrate this in Python using something like scipy.integrate.quad, but this is extremely slow, and this is an integral I need to compute many times over. Does anyone know of a fast way I can approximate an integral of this form?

One idea: if I plot the integrand for some particular values of $q_i$'s, I get something like this: enter image description here Of course, this looks like a Gaussian. But I'm unable to see the reformulation of the integral that would put it in the form of a Gaussian integral. Of course, if this really were a Gaussian integral in disguise, then computing the definite integral is easy with some approximations of the error function.

EDIT: This certainly cannot be reduced to a Gaussian integral. A particular choice of $q_i$'s and $a$ gives the following asymmetric plot of the integrand: enter image description here

  • 1
    $\begingroup$ Are $a$, $b$, $c$, and $d$ the $q_i$? $\endgroup$
    – nicoguaro
    Jun 24, 2021 at 0:08
  • 1
    $\begingroup$ To speed it up I would try a very simple trapezoidal rule (np.trapz) (and compare to quad to check the result / figure out how many points you need). One can speed it up further (if doing many similar integrals) by precomputing the exponential (after a change of variables) at the sample points. This would avoid most of the costly operations. $\endgroup$
    – Winther
    Jun 24, 2021 at 22:47
  • 1
    $\begingroup$ Yeah very often the simplest solution is "best" (well easy and good enough). Just be aware that using it you are giving up the accuracy control (which are built in to quad) and you will be responsible yourself for ensuring the accuracy. So be sure to check the convergence of your results wrt the number of points on some cases before using it. $\endgroup$
    – Winther
    Jun 24, 2021 at 23:27
  • 2
    $\begingroup$ You also have the option of using a fixed Gauss quadrature: scipy.integrate.fixed_quad $\endgroup$
    – nicoguaro
    Jun 25, 2021 at 23:12
  • 1
    $\begingroup$ A fixed Gauss-Hermite quadrature centered at the peak of your function should works. In Python you can find this quadrature in the quadpy library. $\endgroup$ Jun 29, 2021 at 17:19


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.