# Is it possible to make Mathematica reformulate an expression in a more numerically stable way?

I'm writing a numerical optimization, and I'm having a problem with an expression of the form $$e^{-t} (1+\mathrm{erf}(t))$$

The overall shape of the function looks correct, but when $t$ is small, $e^{-t}$ is huge while $(1+\mathrm{erf}(t))$ is very small, and their product is also small. This leads to horrible floating point inaccuracies.

I know, of course, that there are multiple things I can do to remedy this, including scaling my problem so that the values are more reasonably sized. Another is to reformulate the expression to avoid ever computing the huge intermediate values. In the particular example I've given, this is simple: $$e^{-t} (1+\mathrm{erf}(t)) \\ \exp{\left[ \log(e^{-t}) + \log{(1+\mathrm{erf}(t))} \right]} \\ \exp{\left[ -t \log{(1+\mathrm{erf}(t))} \right]}$$ ...which is well behaved for all reasonable values of $t$.

However, in my actual expression, there are various parameters with respect to which I take the derivative. The resulting expressions are hideous and reformulating them by hand is daunting (although tractable).

Is there a way to make Mathematica reformulate an expression while attempting to avoid expressions that will be numerically unstable? I don't expect Mathematica to be automatically aware of which expressions will be problematic, but if I could, for example, simply instruct it to avoid using Exp[] unless it absolutely must, this would be a very useful tool for me (and I suspect for other people working on numerical optimization!).

Note: I am not doing the optimization work in Mathematica. I'm only using Mathematica to help derive analytical gradients for my merit function. Therefore, any features of Mathematica which would eliminate the numerical inaccuracy only in Mathematica doesn't really help me.

-
Eep. Is $t$ positive or negative? Here is a related question. –  Ｊ. Ｍ. May 25 '12 at 14:42
BTW: your function is equivalent to $\exp(-t)\mathrm{erfc}(-t)$; maybe that reformulation might be more useful to you. –  Ｊ. Ｍ. May 25 '12 at 14:47
@J.M.: $t$ can go to negative values, yes. I've actually written my code using $\mathrm{erfc}$, but the values resulting from the $\mathrm{erf}$ term are the same, it doesn't really help. The issue is the huge values from $\mathrm{exp}$. They don't survive to the output anyway, so my goal is to reformulate things such that huge numbers are never generated in the first place. –  Colin K May 25 '12 at 14:56
@J.M.: Matlab however, does have a nice implementation of the scaled complimentary eror function (I'm assuming you mean $e^{t^2} \mathrm{erfc}(t)$) so a reformulation involving this function would be great. I actually hadn't thought of that yet. Thank you! For this particular problem I think I've got enough suggestions to make a solution, but for future use I'd still be interested in how to use MMA for this more effectively. –  Colin K May 25 '12 at 15:25
@Jens: I'm sorry, I should have been more clear about this: The function I've described is a fit function, not the merit function itself. –  Colin K May 25 '12 at 15:43

If you at least know in advance the range in which you will later evaluate, you might consider Taylor series or related approximate forms, using Mathematica to derive such forms and to approximate an error bound. In your example:

In[5]:= func = Exp[-t]*(1 + Erf[t]);
ser = Series[Exp[-t]*(1 + Erf[t]), {t, 0, 4}];
approxpoly = Normal[ser]

Out[7]= 1 + (-1 + 2/Sqrt[Pi])*t + (1/2 - 2/Sqrt[Pi])*t^2 +
(-(1/6) + 1/(3*Sqrt[Pi]))*t^3 + (1/24 + 1/(3*Sqrt[Pi]))*t^4


That is a Taylor polynomial approximation near the origin. We can find the coefficient of the next term to give a first order approximation of the error term.

In[4]:= SeriesCoefficient[Exp[-t] (1 + Erf[t]), {t, 0, 5}] // N
Out[4]= -0.0365428


So the error is around 3/100*t^5 in magnitude near t = zero. This of course assumes the series converges in that region, but in this case at least we know it does (if it did not, your problems would go beyond numerical (in)stability).

At t=1, this error is

In[13]:= func - approxpoly /. t -> 1.

Out[13]= -0.0732347


If you anticipate values of t in that range you might want to use more terms in the polynomial, or else switch before t=1 to a different approximation (this is assuming an error of 7% is more than you want to allow). The table below will give a better idea of the error sizes for t in the range of (-1,1).

In[16]:= Table[func - approxpoly, {t, -1., 1., .1}]

Out[16]= {-0.0239914, -0.0110741, -0.00431147, -0.0012192, \
-0.0000860751, 0.000162876, 0.000119029, 0.0000438515, 7.806*10^-6,
3.05852*10^-7, 0., -4.21962*10^-7, -0.0000151951, -0.000127231, \
-0.000581433, -0.00189768, -0.00499075, -0.0112854, -0.0228166, \
-0.0423096, -0.0732347}


As you would be evaluating outside of Mathematica, you would be substituting the approximations for the function of interest. What I indicate above is an idea of how Mathematica might be utilized to derive and assess them for quality.

-
I would even go further and suggest that OP use PadeApproximant[] to derive approximations, as these can usually give more accurate approximations even slightly far away from the expansion point. –  Ｊ. Ｍ. May 25 '12 at 15:14
For instance, Table[Exp[-t]*(1 + Erf[t]) - g[t], {t, -1., 1., .1}], where g[t_] = PadeApproximant[Exp[-t] Erfc[-t], {t, 0, 3}], gives {0.0175211, 0.0105228, 0.00583002, 0.0029265, 0.00129813, 0.000490888, 0.000150096, 0.0000345791, 6.66391*10^-6, -3.98605*10^-8, 0., 8.10952*10^-9, 7.02389*10^-7, 8.68913*10^-6, 0.0000488437, 0.000178342, 0.000495227, 0.00113726, 0.00227074, 0.00407224, 0.00670781}. –  Ｊ. Ｍ. May 25 '12 at 15:24
@J.M. Thanks. Yes, Pade rational functions are probably better than Taylor polynomials for the stated purpose. –  Daniel Lichtblau May 25 '12 at 15:40

One is reminded of similar problems with $\Gamma$: this function is frequently used in combinatorial expressions in which enormous values almost cancel. Taking logarithms is unstable, but it's possible to compute $\log(\Gamma(z))$ directly with an asymptotic expansion--whence Mathematica offers a separate LogGamma function for this purpose.

Why not use asymptotic expansions for large arguments, then? What is needed is an asymptotic expansion of $\log(\text{erfc}(z))$ as $z \to \infty$ (see the comment by @J.M. and apply $\log$). Mathematica's Limit capability quickly lets us obtain as many terms as we might like. I used it to work out the first few below:

logErfc[z_, n_] :=
With[{t = z^2, c = {1, 5, 37, 353, 4081, 55205, 854197, 14876033}},
-(Log[\[Pi]]/2 + t + (1/2)Log[t]) + Sum[c[[i]] / (i (-2 t)^i), {i, 1, Min[n, Length[c]]}]
]


The n argument lets you choose how many terms to use ($n=0, 1, 2, \dots, 8$). The expression rapidly becomes extremely accurate even for small numbers of terms and modest values of $z$:

ListLogPlot[
Table[Abs[Log[Erfc[N[i, 20]]] - logErfc[i, n]], {n, 2, 8}, {i, 3, 30}],
DataRange -> {3, 30}, PlotRange -> {All, All},
Joined -> True, AxesLabel -> {z, "Absolute log error"}]


To get more coefficients (which you would scarcely need, but this shows how the first eight coefficients were obtained), compute suitable limits of the difference. For instance, to obtain the ninth coefficient, evaluate

Limit[9 (-2 z^2)^9 (Log[Erfc[z]] - logErfc[z, 8]), z -> Infinity]

-
Here is a less complicated way of obtaining a nice asymptotic expansion: Series[Log[Erfc[z]], {z, Infinity, 20}]. One might consider constructing a Padé approximant from the asymptotic series as well, but thanks to a bug within the internal implementation of PadeApproximant[], some trickery is required: Log[1/z] - Log[Pi]/2 - PadeApproximant[Log[1/z] - Log[Pi]/2 - Series[Log[Erfc[z]], {z, Infinity, 10}], {z, Infinity, 6}]. –  Ｊ. Ｍ. May 26 '12 at 2:25
...and if all you want is the c list: CoefficientList[Series[Log[Erfc[Sqrt[z]]], {z, Infinity, 20}] + z - Log[1/Sqrt[z]] + Log[Pi]/2, 1/z] (# 2^#) &@Range[0, 20] –  Ｊ. Ｍ. May 26 '12 at 2:30