Approximation of `erfint`

[MassimoCimmino]

An approximation of the integral of the error function, utilities.erfint(), was proposed by @tblanke in #237.

This issue is to clean up the commits and quantify the impact on accuracy and computational time.

[tblanke]

Hi,

Are there any news regarding this issue?

[MassimoCimmino]

Here are the two implementations I've been comparing:

# Current implementation
def erfint_0(x):
    """
    Integral of the error function.

    Parameters
    ----------
    x : float or array
        Argument.

    Returns
    -------
    float or array
        Integral of the error function.

    """
    return x * erf(x) - 1.0 / np.sqrt(np.pi) * (1.0 - np.exp(-x**2))

# Proposed implementation
def erfint_1(x, tol=1e-8):
    """
    Integral of the error function.

    

    Parameters
    ----------
    x : float or array
        Argument.
    tol : float, optional
        Relative tolerance on the value of erfint.
        Default is 1e-8.

    Returns
    -------
    float or array
        Integral of the error function.

    """
    # Determine the threshold for tolerance
    x_tol = np.maximum(erfinv(1. - tol), np.sqrt(-np.log(tol)))
    # Apply the approximation to all x
    abs_x = np.abs(x)
    y = abs_x - 1. / np.sqrt(np.pi)
    # Evaluate erfint for x < x_tol
    idx = np.less(abs_x, x_tol)
    x_low = x[idx]
    y[idx] = x_low * erf(x_low) - (1.0 - np.exp(-x_low*x_low)) / np.sqrt(np.pi)
    return y

The proposed modified implementation only seems to provide better performance when a sufficient portion of the arguments are above the threshold (x_tol):

All values below x_tol

>>> num = 1000
>>> x_tol = np.maximum(erfinv(1. - 1e-8), np.sqrt(-np.log(1e-8)))
>>> x = np.random.random(num)
>>> %timeit erfint_0(x)
22 µs ± 70.3 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
>>> %timeit erfint_1(x)
35.2 µs ± 68.9 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

25% of values above x_tol

>>> y = x + x_tol - 0.75
>>> %timeit erfint_0(y)
36.4 µs ± 87.4 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
>>> %timeit erfint_1(y)
43.7 µs ± 47 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

50% of values above x_tol

>>> y = x + x_tol - 0.50
>>> %timeit erfint_0(y)
36.3 µs ± 78.4 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
>>> %timeit erfint_1(y)
37.4 µs ± 148 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

75% of values above x_tol

>>> y = x + x_tol - 0.25
>>> %timeit erfint_0(y)
36.3 µs ± 86.2 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
>>> %timeit erfint_1(y)
29.5 µs ± 72.1 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

100% of values above x_tol

>>> y = x + x_tol
>>> %timeit erfint_0(y)
36.2 µs ± 68.9 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
>>> %timeit erfint_1(y)
15.4 µs ± 30.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

It seems that the break even point is around 25% (on my machine). I am hesitant about these changes since some of the cases do run slower. (Updated on 2024-05-23)