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

Question

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.

Answer 1

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.

Answer 2

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]