Efficiently approximating the integral of $\operatorname{erf}(x y) \exp(-x^2)$

Question

I need to efficiently approximate this integral which represents the Gaussian-weighted area of a right triangle whose three points are the origin plus two points that form a vertical "edge", one of which is the edge's nearest point to the origin at $y=0$ and the other, the corner, which is at the same value of $x$. This integral has two inputs, $x$ which is the distance to the edge and $s$ which is the slope from the origin to the corner point. To help with avoiding high degree components the slope needs be no larger than 1 (if the slope is originally larger than 1 the inputs are swapped and the result used to subtract to $\operatorname{erf}(x)\operatorname{erf}(s)$).

$$\int_0^x \frac{\operatorname{erf}(t s) \exp(-t^2)}{2\sqrt{\pi}}\, dt$$

I tried to approximate it in the $x = [-3 , 3], s = [-1 , 1]$ range with a simple 2D Chebyshev approximation, however to achieve my desired maximum absolute error of about $4\cdot10^{-5}$ I need about 31 coefficients in my 2D polynomial. That's a bit much, ideally I would be able to use fewer operations. I have tried to improve the situation by applying operations to the integral before fitting, such as the reciprocal, dividing by $x$ or $s$ or trying to somehow use the limits to gain some sort of advantage, like making a polynomial to interpolate between the limit at $s=0$ (which is $\frac{1-\exp(-x^2)}{2\pi}$) and the integral at $s=1$ (which is $\frac{{\operatorname{erf}(x)}^2}{8}$), or somehow use the limit at $x=\inf$ (which for my purposes is very close to what we have at $x=3$) which is $\frac{\arctan(s)}{2\pi}$, but this didn't help much in achieving more precision with fewer polynomial coefficients. Perhaps there is something clever that I didn't think of or something completely different than a mere polynomial that could be more efficient, perhaps a way to break the problem down into 1D steps.

Integral to approximate, shown with added contours. Yellow is positive, blue negative. Horizontal axis is $x = [-3,3]$, vertical is $s = [-1,1]$.

I'm aware that I could well segment the approximation into a table of small piecewise 2D polynomials, however due to the need to avoid table lookups and dynamic branching I'd rather keep it all in one piece. This integral is used to closely approximate the ideal graphical result of a concave or convex polygon convolved with a 2D Gaussian kernel by calculating the sum of the Gaussian-weighted areas of triangular subdivisions of the polygon, my current implementation works nicely but needs more precision and hopefully can be made more efficient.

Answer 1

Wikipedia has a good discussion on the so-called error function including some approximation formulas available here.

I have employed one of the series approximation in a spreadsheet calculation with good results.

You must, for your problem, combine the above into a good numerical analysis routine to approximate the value of integral.

Good luck (as it is a process of trial and error).

Answer 2

The problem being $s \neq1$, may be, we could try $$\text{erf}(ts)=\frac{2}{\sqrt{\pi }}\sum_{n=0}^\infty (-1)^n \frac{t^{2 n+1}}{(2 n+1) \, n!} s^{2n+1}$$ which would give $$\int_0^x \frac{\operatorname{erf}(t s) \exp(-t^2)}{2\sqrt{\pi}}\, dt=\frac 1{2\pi}\sum_{n=0}^\infty (-1)^n \frac{\Gamma (n+1)-\Gamma \left(n+1,x^2\right)}{(2 n+1) \Gamma (n+1)}s^{2n+1}$$

Computing the partial sums for $s=\frac 12$ and $x=2$ $$\left( \begin{array}{cc} p & \sum_{n=0}^p \\ 0 & 0.078119959 \\ 1 & 0.072095800 \\ 2 & 0.072853673 \\ 3 & 0.072753041 \\ 4 & 0.072765861 \\ 5 & 0.072764343 \\ 6 & 0.072764508 \\ 7 & 0.072764491 \\ 8 & 0.072764493 \end{array} \right)$$

Edit (after @Joako's answer)

Let $t=\frac T s$ and using Mathematica on Wolfram Cloud gives immediately

$$\int_0^x \frac{\operatorname{erf}(t s) \exp(-t^2)}{2\sqrt{\pi}}\, dt=\frac{\tan ^{-1}(s)}{2 \pi }-\text{OwenT}\left[x\sqrt{2} ,s\right]$$ as already given by @Joako.

Expanded as a series around $x=0$ $$\text{OwenT}\left[x\sqrt{2} ,s\right]=\frac{e^{-x^2} }{2 \pi }\sum_{n=0}^\infty \frac{b_n}{n!} x^{2n}$$ $$b_0=\tan ^{-1}(s) \qquad \text{and} \qquad b_{n}=b_{n-1}+(-1)^{n}\frac{ s^{2 n-1}}{2 n-1}$$

Suppose that we just take the expansion to $O(x^6)$ to make $$F(x)=\frac{e^{-x^2}}{2 \pi }\Bigg[\tan ^{-1}(s)+ \left(\tan ^{-1}(s)-s\right)x^2+\frac{s^3-3 s+3 \tan^{-1}(s)}{6} x^4 \Bigg]$$

Computing $$\Phi(s)=\int _{-3}^{+3} \Big[\text{OwenT}\left[x\sqrt{2},s\right]-F(x) \Big]^2\,dx$$

$$\left( \begin{array}{cc} s & \Phi(s) \\ 0.1 & 4.438\times 10^{-19} \\ 0.2 & 6.694\times 10^{-15} \\ 0.3 & 1.709\times 10^{-12} \\ 0.4 & 8.020\times 10^{-11} \\ 0.5 & 1.467\times 10^{-9} \\ 0.6 & 1.469\times 10^{-8} \\ 0.7 & 9.673\times 10^{-8} \\ 0.8 & 4.688\times 10^{-7} \\ 0.9 & 1.800\times 10^{-6} \\ 1.0 & 5.759\times 10^{-6} \end{array} \right)$$

Using the expansion up to $n=7$, the results are

$$\left( \begin{array}{cc} s & \Phi(s) \\ 0.1 & 1.769 \times 10^{-35} \\ 0.2 & 3.047 \times 10^{-30} \\ 0.3 & 2.554 \times 10^{-24} \\ 0.4 & 3.685 \times 10^{-20} \\ 0.5 & 5.672 \times 10^{-17} \\ 0.6 & 2.102 \times 10^{-14} \\ 0.7 & 2.909 \times 10^{-12} \\ 0.8 & 1.957 \times 10^{-10} \\ 0.9 & 7.598 \times 10^{-9} \\ 1.0 & 1.918 \times 10^{-7} \end{array} \right)$$

Update

We can easily build $P_m(x)$, the $[2m,2m]$ Padé approximant of $$2 \pi e^{x^2} \text{OwenT}\left[x\sqrt{2} ,s\right]$$ which is, at least, $O(x^{4m+2})$. Have a look here for an example

Computing again the norm $$\Phi_m(s)=\int _{-3}^{+3} \Big[\text{OwenT}\left[\sqrt{2} x,s\right]-P_m(x) \Big]^2\,dx$$ and reporting the results in a logarithmic scale $$\left( \begin{array}{cccc} s & \log_{10}\Big[\Phi_2(s)\Big]& \log_{10}\Big[\Phi_3(s)\Big]& \log_{10}\Big[\Phi_4(s)\Big] \\ 0.1 & -28.98 & -38.71 & -\infty \\ 0.2 & -22.49 & -29.85 & -36.99 \\ 0.3 & -18.82 & -24.84 & -30.62 \\ 0.4 & -16.34 & -21.46 & -26.29 \\ 0.5 & -14.53 & -18.99 & -23.12 \\ 0.6 & -13.14 & -17.11 & -20.70 \\ 0.7 & -12.04 & -15.63 & -18.79 \\ 0.8 & -11.16 & -14.45 & -17.25 \\ 0.9 & -10.43 & -13.48 & -16.00 \\ 1.0 & -9.820 & -12.69 & -14.95 \end{array} \right)$$ which seems to be much better.

Update We can also write $$\frac{\tan ^{-1}(s)}{2 \pi }-\text{OwenT}\left[x\sqrt{2} ,s\right]=\frac{s x^2}{2 \pi } e^{-x^2}\Bigg[1+3\sum_{n=1}^\infty (-1)^n\frac 1{n!} \frac{a_n}{b_n} x^{2n}\Bigg]$$ where the first $b_n$ are $$\{2,15,140,525,6930,105105\}$$ and the $a_n$ $$\left( \begin{array}{cc} 0 & 1 \\ 1 & s^2-3 \\ 2 & 3 s^4-5 s^2+15 \\ 3 & 15 s^6-21 s^4+35 s^2-105 \\ 4 & 35 s^8-45 s^6+63 s^4-105 s^2+315 \\ 5 & 315 s^{10}-385 s^8+495 s^6-693 s^4+1155 s^2-3465 \\ 6 & 3465 s^{12}-4095 s^{10}+5005 s^8-6435 s^6+9009 s^4-15015 s^2+45045 \end{array} \right)$$

Answer 3

This is just a partial answer:

Here on section $4.3$ eq. $1$ is stated that: $$ \int_{0}^{a} e^{-x^2}\text{erf}(x)\,dx = \frac{\sqrt{\pi}}{4}\Big(\text{erf}(a)\Big)^2$$

So at least for $s=1$ you can use the square of an $\text{erf}$ function (as you mentioned).

Since for other values of $s$ you will have a projection of an always positive Gaussian exponential over different values for $\text{erf}(ts)$ which is a monotonically increasing function, but you won´t be able to solve them in closed form (at least with the mentioned formula), because integration by parts will always left different arguments on the exponential and the error function: $$ \frac{\partial}{\partial x}\Big(\frac{1}{4}\text{erf}(x)\text{erf}(sx) \Big) = \frac{s\,\text{erf}(x)\,e^{-s^2x^2}-e^{-x^2}\,\text{erf}(sx)}{2\sqrt{\pi}}$$ In you shoes I could try with a weighted versions of the same squared error functions or something like $\frac{1}{4}\text{erf}(x)\text{erf}(sx)$ an then try to fit the difference function with a Taylor series or Padé Aproximant of low orders (like this as example)... at least Wolfram-Alpha gives a related-alike answer (but it also give a "$T$" function which I don´t know what it is, but maybe it will be easier than fitting the whole thing).

---

Added later:

Following the comments, since it looks like that the unknown function on Wolfram-Alpha was the Owen's T function, I search about it because of curiosity (I remember see it in some communications papers years ago at university), and by luck I found that surely it is Owen´s T function because of what it is said here on point $3$: $$ \frac{\partial}{\partial z}\Big( T[z;a]\Big) = -\frac{e^{-\frac{z^2}{2}}}{2\sqrt{2\pi}}\text{erf}\left(\frac{az}{\sqrt{2}}\right)$$ $$ T[z;a] = \frac{1}{2\sqrt{2\pi}}\int\limits_{-z}^{+\infty} e^{-\frac{t^2}{2}}\text{erf}\left(\frac{at}{\sqrt{2}}\right)\,dt $$ where it can be easily seen the relation with the integral of the question.

Here some ways to numerically calculate the Owen's T function are shown even with code, and some properties are also listed, like $T[-z;a] = T[z;a]$.

Here I have delete a section that was mistaken.

Applying the fundamental theorem of calculus to the differential formula, and that also $T[0,a] = \frac{\arctan(a)}{2\pi}$ (see here), it looks like a proper development is: $$ \int\limits_{0}^{x} \frac{e^{-t^2}\text{erf}(st)}{2\sqrt{\pi}}\,dt = \frac{\arctan(s)}{2\pi}-T[\sqrt{2}\,x; s]$$ Which looks to work at least for real positive $\mathbf{x,\,s}$ on Wolfram-Alpha, so forget about the previous result (droped since $T[0,a]\neq 1/2 \,\forall a$ - already deleted section), but it still starts integrating from some value different from zero, so keep caution. Also, $T[x,-a] = - T[x,a]$ so the result I give could behave differently from the result of Wolfram-Alpha for $s<0$, so I suggest you to keep using their result.

By the way, on Wolfram-Alpham the Owen´s T function is used as $\text{OwenT[x,a]}$.

I hope at least the founded formulas helps you. Good luck.

---

Last attempt Previously I give an attempt that indeed has a sign mistake that drops the results (as I was afraid). Since some properties was already typed and could help someone to answer the question, I will left them.

Using what is show here: $$T[h;a] + T\left[ah; \frac{1}{a}\right] = \begin{cases} \frac{1}{2} \Phi(h)+\frac{1}{2} \Phi(ah)-\Phi(h)\Phi(ah),\, \text{if}\,s\geq 0 \\ \frac{1}{2} \Phi(h)+\frac{1}{2} \Phi(ah)-\Phi(h)\Phi(ah)-\frac{1}{2} ,\, \text{if}\,s < 0 \end{cases}$$ where the $\Phi(x)$ is the standard Gaussian cumulative distribution $\Phi(x) = \frac{1}{2}\left\{1+\text{erf}\left(\frac{x}{\sqrt{2}}\right) \right\}$,

Noting that now using $h=\sqrt{2}x$ and $a=s$ then: $$\frac{1}{2} \Phi(\sqrt{2}x)+\frac{1}{2} \Phi(\sqrt{2}sx)-\Phi(\sqrt{2}x)\Phi(\sqrt{2}sx) = \frac{1}{4}\left(1-\text{erf}(x)\text{erf}(sx)\right)$$

You will have again that (here was where I made the sign mistake): $$ \begin{array}{r c l} \frac{\partial}{\partial x}\Big(\frac{1}{4}\text{erf}(x)\text{erf}(sx) \Big) & = & \frac{s\,\text{erf}(x)\,e^{-s^2x^2}}{2\sqrt{\pi}}-\frac{\partial}{\partial x}\Big(T[\sqrt{2}x;s]\Big) \\ & = & -\frac{\partial}{\partial x}\Big(T[\sqrt{2}sx;\frac{1}{s}]\Big)-\frac{\partial}{\partial x}\Big(T[\sqrt{2}x;s]\Big) \end{array} $$

So I hardly believe you will escape from the Owen's T functions, at least with a pencil. There are other different expressions for the Owen's T function on this Wiki but I didn´t find how to made simpler expression neither, but maybe they could be easier to become a series, just in case.

As a completely different scope, since Owens' T function is a probability distribution, a things I have used in the past (for other "weird" curves), is to fit another distribution for which the software I use already have calculation packages, one distribution which have result highly flexible is the three parameters Generalized Gamma Distribution, maybe is suitable to approximate the Owen's T functions (but since your requirements of precision are highly accurate, I am no sure if it is going to work for you). Some software like R already have curve fitting for this distribution, maybe it is going to be fast to give it a try.

---

Answering to a comment I made below:

Me: Did you tried series for $T\left[\sqrt{2}sx; \frac{1}{s}\right]$? They left only constants inside the $\text{erf}(x)$ functions when tested on point different from $x=0$, at least they look a bit easier to work with.

OP: I'm not sure what that would develop into.

I am not really sure about this, but it is why I have asked the question. Since you are require an accurate a series representation, that also is efficient in the calculation times, I was thinking in exploring series from points $x \neq 0$.

Since you already have that for positive $x,s$: $$ \int\limits_{0}^{x} \frac{e^{-t^2}\text{erf}(st)}{2\sqrt{\pi}}\,dt = \frac{\tan^{-1}(s)}{2\pi}-T[\sqrt{2}\,x; s]$$ The power series of the integral on the variable $x$ are affine/proportional to the series expansions of $T[\sqrt{2}x;s]$, so as you mentioned knowing it didn´t solve too much for your ambitions.

Playing with Wolfram-Alpha, the $\frac{\tan^{-1}(s)}{2\pi}-T[\sqrt{2}x;s]$ "looks similar" to the $\text{erf}(x)$ from positive $x$: not at zero, but near zero is like a linear function, and far from zero is like a constant, so to fit something in the place where it bends (the harder part), I was thinking on evaluating the series at $x\neq 0$, as example, I have arbitrarily choose $x=1$ (to avoid difficult squared terms and with luck, making easier factorization).

With this, on Wolfram-Alpha the series expansion (Taylor) of $\frac{\tan^{-1}(s)}{2\pi}-T[\sqrt{2}x;s]$ on the variable $x$ at $x=1$ is given by a polynomial $f(s)+\mathbb{P}(x,s,\text{erf(s)})$.

Since the $\text{erf}(x)$ function is known for having slow converging series for $x>1$ (see on Wiki), I was trying somehow to remove part of its weirdness before applying the series expansion.

I am not completely sure about the following formula because of sign issues with the variable $s$ (here a different convention than Wikipedia for the signs of the argument are used I think, this as example), but since it is not the main argument for this partial answer, for now supposed it is true: applying the properties of the Owen's T function $T[h;a]+T[ah;1/a]$ listed above, the integral could be represented as: $$ \int\limits_{0}^{x} \frac{e^{-t^2}\text{erf}(st)}{2\sqrt{\pi}}\,dt = \frac{\tan^{-1}(s)}{2\pi}+\frac{1}{4}\text{erf}(x)\text{erf(sx)}-\frac{\text{sgn}(s)}{4}-T\left[\sqrt{2}sx; \frac{1}{s}\right]$$ Which looks more similar to the answer of Wolfram-Alpha, but they are not identical, so the issues on the sign of $s$ could be messing it up again, so caution if you use it.

Since more weird functions (non-linear for our purpose) have popped out from the Owen's T function (for the formula I give or also for the W-A answer), I was reviewing if this other version of the Owen´s T function was a bit easier for working with it: when viewing the the series expansion (Taylor) of $T\left[\sqrt{2}sx;\frac{1}{s}\right]$ on the variable $x$ at $x=1$ on Wolfram-Alpha, it results it is given by a polynomial $T\left[\sqrt{2}s;\frac{1}{s}\right]+\mathbb{P}(x,s,\text{erf(1)}) \equiv T\left[\sqrt{2}s;\frac{1}{s}\right]+\mathbb{P}(x,s)$ so now the error functions on the coefficients of the series are given only by polynomials of $x$ and $s$, constants $\text{erf}(x=1)$ and a term $T\left[\sqrt{2}s;\frac{1}{s}\right]$ independent of the variable $x$, which maybe could be faster to compute than the series for $T[\sqrt{2}x,s]$.

I don´t know if it really useful, or even if applicable to the Chebyshev polynomials you are using, but since it looks interesting to me I have asked if you have tried it, but more for curiosity than for giving an answer.

Answer 4

I'll keep track of my current best approach. In the question I mentioned my best answer being a simple 31-coefficient polynomial. I just found something a bit better for the same precision, an 18-coefficient polynomial: $$z = (((((-1.6663128{\times}10^{-5}s^2 + 5.104393{\times}10^{-6})x^2 + 0.0005496131s^2 - 5.30433{\times}10^{-5})x^2 + (0.0001584783s^2 - 0.00741157237)s^2 - 0.0018265954)x^2 +(-0.003881739s^2 + 0.0523013844)s^2 + 0.04582956)x^2 + ((-0.00368418s^2 + 0.03692744)s^2 - 0.1996548)s^2 - 0.50360028)x^2 + ((-0.0012717s^2 - 0.0101518)s^2 + 0.0109084)s^2 - 1.836892$$

to which a few operations are applied:

$$\int_0^x \frac{\operatorname{erf}(t s) \exp(-t^2)}{2\sqrt{\pi}}\, dt \approx x^2 s \exp(z)$$

The polynomial was fitted to this function $g(x,s)$:

$$f(x,s) = \int_0^x \frac{\operatorname{erf}(t s) \exp(-t^2)}{2\sqrt{\pi}}\, dt$$ $$g(x,s) = \log(\frac{f(\sqrt x,s)}{xs})$$ in the range $x=[0.4,9.2]$ and $s=[-1,1]$ with 2D Chebyshev polynomials with Chebyshev coefficients under $3{\times}10^{-4}$ culled.

Answer 5

Consider a sum of squared quadratic exponentials approximation to the error function, including linear terms in the quadratics: $$ \widehat{\Phi}_M(t;\,\{a_i,b_i,c_i\})=1-\sum_{i=1}^{M} c_i \exp(-a_i t^2+2b_i t),~t\geq 0 $$ which has $3M$ parameters to be estimated, namely the $\{a_i,b_i,c_i\}$, $i=1,\ldots,M$ with $a_i>0$ $\forall i$. The approximation is valid for nonnegative t, then extended to the real line using the odd symmetry of the error function ${\rm erf}(-t)=-{\rm erf}(t)$. Note that this approximator can be made exact at $t=0$.

Leaving aside the question of the approximation parameters, we can plug in the appoximator into the integral in question $$ I(x,s)=\frac{1}{2\sqrt{\pi}}\int_0^x {{\rm erf}}(t s) \exp(-t^2)\, dt $$ to obtain $$ I(x,s)\approx \frac{1}{2\sqrt{\pi}}\int_0^x\exp(-t^2)\, dt-\frac{1}{2\sqrt{\pi}} \sum_{i=1}^{M} c_i \int_0^x \exp(-a_is^2t^2+2b_ist)\exp(-t^2)\, dt $$ which simplifies to $$ I(x,s)\approx \frac{1}{4}{\rm erf}(x) -\frac{1}{2\sqrt{\pi}} \sum_{i=1}^{M} c_i \int_0^x \exp(-\alpha_i^2t^2+2\beta_it)\, dt $$ where we have defined $\alpha_i^2=1+a_i s^2$ and $\beta_i=b_i s$. The integral above is standard, even in the indefinite case, and according to 2.33(1) of Gradshteyn and Ryzhik's "Table of Integrals, Series and Products," we have (for $\alpha\neq 0$, with $\alpha>0$ here): $$ \int \exp(-\alpha^2t^2+2\beta t)\, dt=\frac{\sqrt{\pi}}{2\alpha}\exp\left(\frac{\beta^2}{\alpha^2}\right){\rm erf}\left(\alpha t-\frac{\beta}{\alpha}\right) $$ Therefore we get an approximator of the form $$ I(x,s)\approx \frac{1}{4}{\rm erf}(x) -\sum_{i=1}^{M} \frac{c_i}{4\alpha_i} \exp\left((\beta_i/\alpha_i)^2\right) \left[ {\rm erf}\left(\alpha_ix-\frac{\beta_i}{\alpha_i}\right)+ {\rm erf}\left(\frac{\beta_i}{\alpha_i}\right) \right] $$ This approximator only requires $M$ exponential and $2M+1$ error function evaluations and is a direct approximation - not piecewise for the entire parameter region in $x$ and $s$.

Now the interesting part: how to obtain the coefficients $\{a_i,b_i,c_i\}$, $i=1,\ldots,M$? This is the subject of a current paper I have in review and I can't give full details here at the moment, suffice to say it is not difficult - a simple gradient descent is used to optimise them. I have obtained the coefficients for the $M=4$ case and the fit to the error function over the range $0\leq x\leq 4$ has a peak absolute relative error (ARE) of about $7E-4$ for $x\in[0,4]$, which immediately extends to $[-4,4]$ due to erf being an odd function. This should give a better accuracy under integration, maybe around $1E-5$. The $M=4$ coefficients are listed below. If these do not provide the required accuracy for the integral, one could move to $M=5$, which would probably do the trick.

$a_1,\ldots,a_4$: $+1.102149$, $+0.602149$, $+0.802149$, $+0.302149$

$b_1,\ldots,b_4$: $-0.738479$, $-0.738479$, $-0.638479$, $-0.238479$

$c_1,\ldots,c_4$: $-0.656344$, $-8.65439E\!-\!2$, $+1.742885$, $2.31093E\!-\!6$

The context of this work is the study of low-dimensional generative adversarial networks (GANs) based on error function integrals. Please refer to https://www.researchgate.net/publication/356815736_Convergence_and_Optimality_Analysis_of_Low-Dimensional_Generative_Adversarial_Networks_using_Error_Function_Integrals