5
$\begingroup$

I have a function $g(E)$ that is defined by a very complicated expression but that only involves built-in functions and integrations. I would like to define it in the absolute most efficient way possible. First, I will define this function mathematically, and then I will show two ways I have implemented this.


Function Definition.

My function is actually multi-variable, i.e. it's really $g(E; x,p_z,a)$, and defined as follows:

$$\begin{align} \beta(x)&\equiv 1-x \tag{1}\\ p_1^2&\equiv\left(k_4+\frac{1}{2}E\right)^2+k_{\perp}^2+(xp_z)^2\tag{2}\\ p_2^2&\equiv\left(k_4-\frac{1}{2}E\right)^2+k_{\perp}^2+(\beta(x)p_z)^2\tag{3}\\ M(k)&\equiv\frac{M_0}{(1+k^2)^{a}}\tag{4}\\ E^2_{\pi}&\equiv p_z^2+m^2\tag{5}\\ & \end{align}$$

$$\begin{align} g(E;x,p_z,a)&=\frac{1}{E^2+E_{\pi}^2}\times\\ &\int_{-\infty}^{\infty}dk_4\int_0^{\infty}d(k_{\perp}^2)\frac{\sqrt{M(p_1)M(p_2)}\left(\beta(x)M(p_1)+xM(p_2)\right)}{(p_1^2+M(p_1)^2)(p_2^2+M(p_2)^2)} \end{align}$$

where I have suppressed much functional dependence notation for brevity. I am interested in the following ranges of parameters: $$\begin{align} -1<&x<3\tag{7}\\ 0<&p_z<\infty\\ 1<&\alpha<3 \end{align}$$


Implementations.

Here is my first "natural" implementation where I simply make function definitions as I did just above.

(*Implementation #1: Slow*)

ClearAll[B,M0,p1,p2,M,N1,N2,D1,fint,Epi,g];

M0 = 1.; (*set this equal to 1 for simplicity*)
B[x_?NumericQ] := 1. - x;
p1[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
   pz_?NumericQ] := (k4 + 0.5*En)^2 + kp + (x*pz)^2;
p2[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
  pz_?NumericQ] := (k4 - 0.5*En)^2 + kp + (B[x]*pz)^2;
M[k_?NumericQ, a_?NumericQ] := M0/(1. + k)^a;

N1[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
   pz_?NumericQ, a_?NumericQ] := 
  Sqrt[M[p1[k4, kp, x, En, pz], a]*
    M[p2[k4, kp, x, En, pz], a]]; (*numerator term #1*)

N2[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
    pz_?NumericQ, a_?NumericQ] := B[x]*
     M[p1[k4, kp, x, En, pz], a] + 
    x*M[p2[k4, kp, x, En, pz], a]; (*numerator term #2*);

D1[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
  pz_?NumericQ, a_?
   NumericQ] := (p1[k4, kp, x, En, pz] + 
    M[p1[k4, kp, x, En, pz], a]^2)*(p2[k4, kp, x, En, pz] + 
    M[p2[k4, kp, x, En, pz], a]^2); (*denominator term #1*)

fint[k4_?NumericQ, kp_?NumericQ, x_?NumericQ, En_?NumericQ, 
   pz_?NumericQ, a_?NumericQ] := (
  N1[k4, kp, x, En, pz, a]*N2[k4, kp, x, En, pz, a])/
  D1[k4, kp, x, En, pz, a]; 

Epi[pz_?NumericQ] := Sqrt[pz^2 + 0.24^2];

g[x_?NumericQ, En_?NumericQ, pz_?NumericQ, a_?NumericQ] := 
  (1./(En^2+Epi[pz]^2))*NIntegrate[
   fint[k4, kp, x, En, 
    pz, a], {k4, -\[Infinity], \[Infinity]}, {kp, 
    0, \[Infinity]}, 
   Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 15000}];

One function call on my computer takes approximately $0.3$ seconds.

g[0.5, 0, 3, 2] // AbsoluteTiming
{0.313269, 0.000160885}

Here is my second (faster) implementation where I forgo all those sub-function definitions and simply make one humongous definition for $g(E)$.

(*Implementation #2: Faster*)

ClearAll[Subscript[M, 0], g];
Subscript[M, 0] = 1.;

g[x_?NumericQ, En_?NumericQ, pz_?NumericQ, a_?NumericQ] := 
  1./(En^2 + pz^2 + 0.24^2)*
   NIntegrate[((Sqrt[
       Subscript[M, 
        0]/(1. + (k4 + 0.5*En)^2 + kp + (x*pz)^2)^a*Subscript[
        M, 0]/(1. + (k4 - 0.5*En)^2 + 
          kp + ((1. - x)*pz)^2)^a])*((1. - x)*Subscript[M, 
        0]/(1. + (k4 + 0.5*En)^2 + kp + (x*pz)^2)^a + 
       x*Subscript[M, 
        0]/(1. + (k4 - 0.5*En)^2 + 
          kp + ((1. - x)*pz)^2)^a))/(((k4 + 0.5*En)^2 + 
       kp + (x*pz)^2 + (Subscript[M, 
        0]/(1. + (k4 + 0.5*En)^2 + kp + (x*pz)^2)^a)^2)*((k4 - 
         0.5*En)^2 + 
       kp + ((1. - x)*pz)^2 + (Subscript[M, 
        0]/(1. + (k4 - 0.5*En)^2 + 
          kp + ((1. - x)*
            pz)^2)^a)^2)), {k4, -\[Infinity], \[Infinity]}, \
{kp, 0, \[Infinity]}, 
    Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 15000}];

Now the same single function evaluation is about 22 times faster.

g[0.5, 0, 3, 2] // AbsoluteTiming
{0.014624, 0.000160885}

However, I need this to be faster. I would like to implement this function in literally the most efficient way possible. How can I make my implementation of $g(E)$ more efficient?


Ultimate Goal.

My end goal is to evaluate the Fourier transform of this function $g(t)$, which requires numerous function calls/calculations, and study the large $t$ limit as a function of $x$ and $p_z$.

$$I(t;x,p_z,a)=\int_{-\infty}^{\infty}dE\,\cos(Et)g(E;x,p_z,a) $$

where I have only written $\cos(Et)$ instead of $e^{it}$ because I already know $g(E)$ is an even function. Right now a single Fourier transform evaluation takes approximately $40$ seconds.

$\endgroup$
  • $\begingroup$ The remaining time is likely not due to a bad definition. It's for computing the (probably singular) integral. That needs time. Certainly, it can be sped up but that would need quite a lot of hand tuning. $\endgroup$ – Henrik Schumacher Jun 25 '18 at 8:49
  • $\begingroup$ Your first attempt is slow because you use delayed evaluation (:=, i.e. SetDelayed and NumericQ), which means it has to evaluate all of those right hand sides every time, once it gets to specific numbers. If you remove those except for the ones in the definition of g then you'll see the speedup but the code will remain readable. $\endgroup$ – KraZug Jun 25 '18 at 12:29
  • 1
    $\begingroup$ Playing around with the method options, "LocalAdaptive" gives the same result in about 50-75% of the time for some sample values. If you also lower the AccuracyGoal you can speed it up another 10-fold, the obvious trade-off. $\endgroup$ – KraZug Jun 25 '18 at 12:48

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.