7
$\begingroup$

Let us say that I have a function like so:

def f(x):
    return g(x)*h(x)

Now, g(x) and h(x) are such that when g(x) produces a result large enough to cause overflow errors, h(x) produces a number that is approximately zero. So, the result provided by f(x) should be $0$.

Would it be enough to write f(x) as follows:

def f(x):
    h_x = h(x)
    if abs(h_x) < 1e-9:
        return 0
    else:
        return g(x)*h_x

Or, is there a neater, more general way to write such a function? If one wants an example of functions f, g and h that behave as described, consider:

\begin{align} h(x) &= \textrm{erf}(a) - \textrm{erf}(x) \\ g(x) &= e^{x^2} \\ f(x) &= g(x)h(x) \\ x &= 30 \\ a &= 30 \end{align}

$\endgroup$
  • 8
    $\begingroup$ It's unclear why you think the product $g(x)*h(x)$ should be approximately zero, just because $h(x)$ is close to zero. After all you open by saying that in these cases $g(x)$ is very large (causing overflow). Perhaps the crux of the matter is how to compute $g(x)$ to avoid overflow. Can we instead compute $\ln g(x)$? Then the rest of the computation could be $f(x) = exp(\ln g(x) + \ln h(x))$, and this would likely be more numerically stable. $\endgroup$ – hardmath Jul 8 '15 at 21:58
  • $\begingroup$ How do you want this to work for $g(x)=e^{x^2}$, $h(x) = e^{-x^2}$? $\endgroup$ – Kirill Jul 8 '15 at 22:22
  • $\begingroup$ @Kirill Good question. I would want it to return 1 in that case, but my suggested f(x) would not do that. $\endgroup$ – user89 Jul 8 '15 at 22:30
  • 1
    $\begingroup$ Sure, if $h(x)$ is precisely zero (or if $g(x)$ is), then their product is precisely zero. Otherwise the relative accuracy of $g(x)\cdot h(x)$ depends on getting relative accuracy in both factors. Aiming at the computation of the logarithms, rather than $g(x),h(x)$ themselves, is a fairly common way of managing the extremes of magnitudes. I would adopt just such a strategy if $h(x) = \operatorname{erf}(a) - \operatorname{erf}(x)$ for large $x \gt a$. $\endgroup$ – hardmath Jul 9 '15 at 0:03
  • 1
    $\begingroup$ I'd be glad to. For the specific case where you want $h(x) = \operatorname{erf}(a) - \operatorname{erf}(x)$, for large $a,x$, a possibility is to compute instead $h(x) = \operatorname{erfc}(x) - \operatorname{erfc}(a)$. Using the complementary error function permits us to take the difference of two very small numbers rather than than two numbers both very close to $1$. I'll discuss this example at greater length in an answer. $\endgroup$ – hardmath Oct 2 '15 at 12:17
6
+50
$\begingroup$

Initially the Question focused on a product computation, f(x)*g(x), where one factor is very small and the other is very large (perhaps large enough to cause overflow in floating point representation).

I commented that it might avoid overflow (resp. underflow) if instead we computed (numerical approximations to) $\ln f(x)$ and $\ln g(x)$ and added them. In exact computation, assuming both $f(x), g(x) \gt 0$:

$$ f(x) g(x) = \exp ( \ln f(x) + \ln g(x) ) $$

In floating point computation the relative precision of f(x)*g(x) would then correspond to the absolute precision of log(f(x)) + log(g(x)). In this sense we have only moved the problem around. As the OP noted in Comments, the numerical evaluation of $f(x)$ might be exactly zero, despite the true value being only a very small non-zero number.

For example, the OP notes, we might have $f(x) = \operatorname{erf}(a) - \operatorname{erf}(x)$ when $x,a$ are both large. In this case both error function values rapidly approach $1$ from below. A reasonable library implementation will return exactly $1$ for all sufficiently large calling arguments. Thus the difference $f(x)$ will be naively computed as exactly zero, a clear illustration of catastrophic cancellation.

Error function from Wikipedia

Fig. 1: Graph of error function, from Wikimedia

In this case the root cause of loss of significance is well understood, and in coding we can arrange a less naive evaluation of $f(x)$, one that allows us to accurately evaluate $\ln f(x)$.

First note the definition of the complementary error function:

$$ \operatorname{erfc}(x) = 1 - \operatorname{erf}(x) $$

For large $x$ we have $\operatorname{erfc}(x)$ close to zero (instead of $\operatorname{erf}(x)$ being close to one). Consequently:

$$ f(x) = \operatorname{erfc}(x) - \operatorname{erfc}(a) $$

avoids catastrophic cancellation (unless $x,a$ are especially close).

As coders we do have responsibility for keeping track of the sign of $f(x)$, since we can only take a logarithm of a positive real number. This task is easy here, since $\operatorname{erf}(x)$ is monotone increasing, so $f(x) \gt 0$ if and only if $x \lt a$. Should $x = a$, then $f(x) = 0$ and the product $f(x) g(x) = 0$ as well. If $x \gt a$, then $f(x) \lt 0$, but we can account for this and work with $\ln |f(x)|$ to compute $\ln |f(x) g(x)|$, fixing up the sign of $f(x) g(x)$ accordingly.

As a further precaution we will likely want to guard against underflow in evaluating $\operatorname{erfc}(x)$. The exponent range for double precision normal values goes down to binary (power of two) $-1022$. Asymptotically:

$$ \operatorname{erfc}(x) \approx O(1/x) e^{-x^2} $$

Consequently we would be at the brink of exponent underflow when $x \approx 25$.

To avoid this some mathematical libraries (MatLab,SciPy) implement the scaled complementary error function:

$$ \operatorname{erfcx}(x) = e^{x^2} \operatorname{erfc}(x) $$

If this special function is not available in your library of choice, some self-implementation ideas are kicked around here, especially Steven G. Johnson's Faddeeva Package implementation notes and its algorithms references.

Successful handling of this special case (difference of two error function values) does not lead to immediate resolution of all problems with numerical accuracy! It merely demonstrates that when a cause of catastrophic cancellation is foreseen, we have an opportunity to "code around" it. Sadly it is often quite time consuming to work out the alternative evaluation and (especially) to test the various paths of computation to ensure consistency.

$\endgroup$
  • $\begingroup$ I thought you might be interested to know that I am having some trouble understanding how $\textrm{erfcx}$ would help me avoid roundoff errors, so I made a question about it here: scicomp.stackexchange.com/questions/21016/… $\endgroup$ – user89 Oct 15 '15 at 21:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.