# How to address the numerical instability of a function involving erf and an exponential

We need to compute this function:

$$f(x) = x + \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- ( \frac {x} {\sqrt 2 \cdot \sigma_x} )^2}} {1 + erf(\frac {x} {\sqrt 2 \cdot \sigma_x})}$$

We noticed that, for negative values of $$\frac {x} {\sqrt 2 \cdot \sigma_x}$$ lower than a certain value (and different depending on which software we use to compute it, say about $$-5$$ for instance), there is numerical instability (oscillations).

I saw this post suggesting how to achieve better stability, although on a rather different function:

Numerical Instability with this cumulate

Do you think a similar method would work here? And if so, should I express the $$erf$$ as an integral?
Or would you suggest some other method? Taylor expansion? Anything else...?

Thanks

EDIT after evaluation of the responses

Thank you all very much for your feedback!

I think there is a complication here that I had not explicitly foreseen.

First I will make the following substitution, so your alternative expressions of $$erf$$ can be used directly:

$$x = y \cdot \sqrt 2 \cdot \sigma_x$$

$$f(y) = \sqrt 2 \cdot \sigma_x \cdot y + \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)}$$

Now, the theory we used to derive this equation requires that:

$$f(y) > 0, \forall y$$
$$\lim_{y \to -\infty} f(y) = 0$$

So, we expect $$f(y)$$ to be very slowly approaching $$0$$ as $$y$$ becomes more negative; but $$f(y)$$ should never be negative.

In fact, if I plot the second term of $$f(y)$$, it looks like, as $$y$$ becomes more and more negative:

$$\sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)} \approx - \sqrt 2 \cdot \sigma_x \cdot y$$

This is in substantial agreement with Claude's approximation for the $$erf$$ part.
I.e., using:

$$\frac{e^{-y^2}}{1+\text{erf}(y)} \approx -\sqrt{\pi } \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg]$$

I get:

$$\sqrt {\frac 2 \pi} \cdot \sigma_x \cdot \frac {e^{- y^2}} {1 + erf(y)} \approx \sqrt {\frac 2 \pi} \cdot \sigma_x \cdot (-\sqrt{\pi }) \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg] =$$
$$= - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg]$$

which indeed approaches $$- \sqrt 2 \cdot \sigma_x \cdot y$$ when $$y$$ becomes very negative.

With the above approximation:

$$f(y) \approx \sqrt 2 \cdot \sigma_x \cdot y - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[1+\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg] =$$ $$= - \sqrt {2} \cdot \sigma_x \cdot \,y\Bigg[\frac{1}{2 y^2}-\frac{1}{2 y^4}+\frac{5}{4 y^6}+O\left(\frac{1}{y^8}\right)\Bigg]$$

Numerically, this 'does the trick', i.e. for $$y \le -4$$ it seems really very close to the function I need to compute, as indicated by Claude. And in fact with the rational function version even up to $$y \le -1$$.

This seems to imply that only polynomial approximations are viable, not the $$erfc$$-based ones, or am I wrong?

Any thoughts?

I would suggest a series expansion for $$x \leq- 4$$ $$f(x)=\frac{e^{-x^2}}{1+\text{erf}(x)}=-\sqrt{\pi } \,x\Bigg[1+\frac{1}{2 x^2}-\frac{1}{2 x^4}+\frac{5}{4 x^6}+O\left(\frac{1}{x^8}\right)\Bigg]$$ For $$x=-5$$, the "exact" value is $$9.03318$$ while this truncated expansion gives $$9.03313$$.

Much better would be the Padé approximants $$f(x)=-\sqrt{\pi } \,x \,\frac{8 x^6+84 x^4+210 x^2+105 } {8 x^6+80 x^4+174 x^2+48 }$$ which, for $$x=5$$, leads to an error of $$2.2\times 10^{-8}$$.

• Thank you Claude, your solution seems to work very well! I am probably going to accept it as the answer. Could you please just comment briefly on how you obtained that expansion? – user6376297 Mar 10 at 13:16
• @user6376297. The first one is obtained by composition of the Taylor series for large negative values of $x$. Taking more terms, it is easy to build the corresponding Padé approximant. The one I produced is equivalent to a Taylor series to $O\left(\frac{1}{x^{13}}\right)$. If you want better, let me know since it just need a couple of minutes to build. Cheers :-) – Claude Leibovici Mar 11 at 2:21
• Thanks, no, I do not need a more precise approximation, I was just trying to understand how you made yours, so for similar cases I could do it myself. Unfortunately I cannot see how that is a Taylor series, as per definition it should only have positive powers of x. en.wikipedia.org/wiki/Taylor_series . Also, I would not know on which point to 'center' the series, as here we say 'very negative x', which for me means make x go to $-\infty$, but then I cannot write Taylor terms with $(x+\infty)^n$... So I am a bit stuck. I will do more research. And I will accept your answer. – user6376297 Mar 11 at 10:05
• OK I found two posts that maybe clarify this. math.stackexchange.com/questions/595426/… and math.stackexchange.com/questions/51770/… . Perhaps there are some substitution and back-substitution steps being taken, not simply a Taylor expansion. The other issue is that the CAS I am using to do this gets the limit wrong for $f(x)$ and its derivatives at $x \to - \infty$. But OK, that can be addressed somehow. – user6376297 Mar 11 at 10:14
• No, tried it out and still does not work for me :( I might post another question about that, unless you can comment on the method you used. – user6376297 Mar 11 at 10:45

In the negatives, $$\text{erf(x)}$$ very quickly tends to $$-1$$ and it is no surprise that you soon reach the limits of the floating-point representation.

You will regain full accuracy by computing

$$\text{erf}(x)+1$$ directly. It is possible that you will obtain better results with $$\text{erfc}(-x)$$, but that depends on the numerical library you use.

Otherwise, you could use the fourth approximation in https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions, and you will see a nice simplification with your numerator.

$$1 + erf(x)$$ suffers from catastrophic cancellation in the negative half plane. This is avoided by replacing it with $$erfc(-x)$$, that is, use of the complementary error function. Most standard math libraries, for example the C++ standard math library, include this as a function erfc.

Many special function libraries also offer the exponentially scaled complementary error function, $$e^{x^{2}} erfc(x)$$, often as a function called erfcx. With that, one can compute $$\frac{e^{-x^{2}}}{1+erf(x)}$$ accurately as 1.0 / erfcx (-x). I also provided a double-precision implementation of erfcx on Stackoverflow.

As I became aware belatedly, there is a second instance of catastrophic cancellation in $$f(x)$$ for negative $$x$$ large in magnitude: $$x$$ and the rest of the expression are similar in magnitude but of opposite sign. Around $$-10^{7}$$ complete loss of accuracy occurs.