How to avoid overflow error in program that computes product of two numbers, such that when one is big enough to cause overflow, other is $0$?

Question

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}

Answer 1

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.

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.