More precision in numpy arrays

Question

I have the following simple function to evaluate.

def f0(wt):
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    return 0.5 * (term1 + term2 + term3)

For small values of wt (order of 1e-4 and below), I seem to have numerical problems in the evaluation of the function. Indeed, the term1 and term3 have very large and almost opposite values, but term2 is very small.

I think I improved things slightly by splitting the sum of the 3 terms into two parts, as showed here

def f1(wt):
    # Split the calculation to have more stability hopefully
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    partial = term1 + term3
    return 0.5 * (partial + term2)

However, for very small but positive values of wt, I think there are still numerical problems. I expect this function to be smooth for any positive value of wt, but, as you can see from the plot attached, at values below 1e-3, there are wild artifacts.

My question is: how can I improve the numerical precision of Numpy, if I am already using the data type float64?

Note: I am on a Windows 10 machine with 64 bits. I have read on other Stack Overflow threads that the class np.float128 is not available.

Full code snippet

import numpy as np
import matplotlib.pyplot as plt

wt = np.logspace(-6, 1, 1000)

def f0(wt):
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    return 0.5 * (term1 + term2 + term3)

def f1(wt):
     # Split the calculation to have more stability hopefully
     term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
     term2 = np.sin(wt)**2
     term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
     partial = term1 + term3
     return 0.5 * (partial + term2)

plt.figure()
plt.loglog(wt, f0(wt), label='f0')
plt.loglog(wt, f1(wt), label='f1')
plt.grid()
plt.legend()
plt.xlabel('wt')
plt.show()

I think the question is not that of improving the precision of Numpy (which ultimately depends on the precision of your machine and underlying math routines as well), but instead focusing on the implementation of this function, which may need annoying special cases e.g. for small values. — Iguananaut
– Iguananaut, Commented Aug 26, 2021 at 15:58

Lukas S · Accepted Answer · 2021-08-27 00:53:13Z

How about you replace the sin and cosin with the first few terms of their Taylor series. Then sympy is able to give you a simple result that is hopefully better suited numerically.

First I slightly change your function so it gives me a sympy expression.

from sympy import *
t = symbols('t')

def f0(wt):
    term1 = (1 + sympy.cos(wt)**2) * (sympy.Rational(1,3) - 2 / (wt)**2)
    term2 = sympy.sin(wt)**2
    term3 = 4 / (wt)**3 * sympy.cos(wt) * sympy.sin(wt)
    return sympy.Rational(1,2)*(term1 + term2 + term3)
expr = f0(t)
expr

Now I replace sin and cos with their taylor polynomials.

def taylor(f, n):
    return sum(t**i/factorial(i) * f(t).diff(t, i).subs(t,0) for i in range(n))

tsin = taylor(sin, 7)
tcos = taylor(cos, 7)

expr2 = simplify(expr.subs(sin(t),tsin).subs(cos(t),tcos))
f1 = lambdify(t, expr2, 'numpy')
expr2

And finally I plot it using exactly your code. Notice that I am using sympys option to make a numpy ufunc.

wt = np.logspace(-6, 1, 1000)
plt.figure()
plt.loglog(wt, f0(wt), label='f0')
plt.loglog(wt, f1(wt), label='f1')
plt.grid()
plt.legend()
plt.xlabel('wt')
plt.show()

Obviously this function is only good around zero and for values between 1 and 10 you should take the original function. But in case you need convincing and don't care that the function with replaced taylor polynomial looks nice you can crank the degree up to 25 making it visually agree with your function at least up until 10.

And you can combine the functions so it calculates the values around zero with my function and the other with yours like this.

def f2(wt):
    cond = np.abs(wt) > 1/10
    return np.piecewise(wt, [cond, ~cond], [f0,f1])

Jérôme Richard · Accepted Answer · 2021-08-26 23:46:45Z

The problem you are facing is catastrophic cancellation and it must not be solved using higher precision as doing so will generally postpone the actual problem. The root of the problem which is a numerical instability must be solved by reformulating the mathematical expression. Note that f1 is a bit better than f0 but the cancellation issue lies in term1 + term3.

By transforming the expression simple development/factorization operations and using trigonometric identities one can get the following function:

def f2(wt):
     sw = np.sin(wt)
     sw2 = np.sin(2*wt)
     return (sw/wt)**2 + 1/3 + (sw2 / wt - 2) / wt**2 + sw**2 / 3

This function is a bit more accurate but still contains a cancellation causing the same issue. This happens because of the expression E = (sw2 / wt - 2) / wt**2 which is the root of the problem. Indeed, np.sin(2*wt) tends towards 2 when wt is near 0. Thus sw2 / wt - 2 is close to 0 and the expression E is numerically unstable because of a close-to-zero value divided by another close-to-zero value. If one can reformulate analytically E to remove the singularity, then the resulting expression will likely be numerically stable. For more information you can look at the sinc function and how to compute an approximation of this function (also available in Numpy).

One simple way to solve this is to use numerical tools like Taylor series. Taylor series can approximate the expression of E close to zero accurately (because of its derivatives). Actually, one can use Taylor series to compute the whole expression and not only E. However, using Taylor series for values close to 1 give inaccurate results. In fact, the accuracy of the method drops very quickly above 1. One solution is to only use the Taylor series for small values.

Here is the resulting implementation:

def f3(wt):
     sw = np.sin(wt)
     sw2 = np.sin(2*wt)
     reference = (sw/wt)**2 + 1/3 + (sw2 / wt - 2) / wt**2 + sw**2 / 3
     # O(13) Taylor series computation used only for near-zero values
     taylor = (  (  4. /        15.) * wt**2
               - ( 29. /       315.) * wt**4
               + ( 37. /      2835.) * wt**6
               - (151. /    155925.) * wt**8
               + (268. /   6081075.) * wt**10
               - (866. / 638512875.) * wt**12)
     # Select the best implementation
     return np.where(np.logical_and(wt >= -0.2, wt <= 0.2), taylor, reference)

This implementation appear to be very accurate in practice (>=12 digits of precision) while being still relatively fast. Here is the result:

Collectives™ on Stack Overflow

More precision in numpy arrays

Full code snippet

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

Full code snippet

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related