Is there a faster way of combining functions containing mpf (mpmath floats) with numpy arrays?

Question

I was having problems with the accuracy of floats in Python. I need high accuracy because I want to use explicitly written spherical bessel functions J_n (x), which deviate (especially for n>5) from their theoretical values at low x values if numpy floats are used (15 precise digits).

I have tried many options, especially from mpmath and sympy, in order to keep more precise numbers. I had problems when combining the accuracy of mpmath inside the functions with numpy arrays, until I knew there was the function numpy.vectorize. Finally I got this solution to my initial problem:

import time
% matplotlib qt
import scipy
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from sympy import *
from mpmath import *
mp.dps=100


#explicit inaccurate
def bessel6_expi(z):
    return -((z**6-210*z**4+4725*z**2-10395)*np.sin(z)+(21*z**5-1260*z**3+10395*z)*np.cos(z))/z**7

#explicit inaccurate 1, computation time increases, a bit less inaccuracy
def bessel6_exp1(z):
    def bv(z):
        return -((z**6-210*z**4+4725*z**2-10395)*mp.sin(z)+(21*z**5-1260*z**3+10395*z)*mp.cos(z))/z**7
    bvec=np.vectorize(bv)
    return bvec(z)

#explicit accurate 2, computation time increases markedly, accurate
def bessel6_exp2(z):
    def bv(z):
        return -((mpf(z)**mpf(6)-mpf(210)*mpf(z)**mpf(4)+mpf(4725)*mpf(z)**mpf(2)-mpf(10395))*mp.sin(mpf(z))+(mpf(21)*mpf(z)**mpf(5)-mpf(1260)*mpf(z)**mpf(3)+mpf(10395)*mpf(z))*mp.cos(mpf(z)))/mpf(z)**mpf(7)
    bvec=np.vectorize(bv)
    return bvec(z)

#explicit accurate 3, computation time increases markedly, accurate
def bessel6_exp3(z):
    def bv(z):
        return -((mpf(z)**6-210*mpf(z)**4+4725*mpf(z)**2-10395)*mp.sin(mpf(z))+(21*mpf(z)**5-1260*mpf(z)**3+10395*mpf(z))*mp.cos(mpf(z)))/mpf(z)**7
    bvec=np.vectorize(bv)
    return bvec(z)

#implemented in scipy, accurate, fast
def bessel6_imp(z):
    def bv(z):
        return scipy.special.sph_jn(6,(z))[0][6] 
    bvec=np.vectorize(bv)
    return bvec(z)

a=np.arange(0.0001,17,0.0001)

plt.figure()
start = time.time()
plt.plot(a,bessel6_expi(a),'b',lw=1,label='expi')
end = time.time()
print(end - start)

start = time.time()
plt.plot(a,bessel6_exp1(a),'m',lw=1,label='exp1')
end = time.time()
print(end - start)

start = time.time()
plt.plot(a,bessel6_exp2(a),'c',lw=3,label='exp2')
end = time.time()
print(end - start)

start = time.time()
plt.plot(a,bessel6_exp2(a),'y',lw=5,linestyle='--',label='exp3')
end = time.time()
print(end - start)

start = time.time()
plt.plot(a,bessel6_imp(a),'r',lw=1,label='imp')
end = time.time()
print(end - start)

plt.ylim(-0.5/10**7,2.5/10**7)
plt.xlim(0,2.0)
plt.legend()

plt.show()

The problem I have now is that just for plotting the explicit, accurate ones, it takes quite a long time (about 31 times slower than the scipy function for mp.dps=100). Smaller dps do not make these processes much faster, even with mp.dps=15, they are still 26 times slower. Is there a way to make this faster?

If you are going to use non-primitive numeric types in `numpy`, you are going to lose your speed gains. Note, from the [docs of `numpy.vectorize`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html) "The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop." — juanpa.arrivillaga, Feb 14 '18 at 18:27
And is there a way not to use these `mpf`s and still have a decent precision? Or (I suppose this may be almost impossible in my case) avoiding using numpy? Thanks btw. — Edi, Feb 14 '18 at 18:39
Yeah, I don't know exactly what would count as "decent precision", but if `np.float64` doesn't work for you, then I think you'd have to use something like this or `decimal.Decimal` — juanpa.arrivillaga, Feb 14 '18 at 18:47
My `numpy` has a `np.float96`, and some may have a `128`. `object` dtype arrays can hold non-numpy objects, but you loose most of the fast compiled functionality. — hpaulj, Feb 14 '18 at 19:00

Paul Panzer · Accepted Answer · 2018-02-15T04:19:50.930

Note that the loss of accuracy you observe near zero comes from the fact that you are subtracting two nearly equal terms both of the form 10395 z^-6 + O(z^-4). As the true value is 1/135135 z^6 + O(z^8) you will lose a factor of ~1.4 x 10^9 z^-12 in accuracy. So if you want to calculate the value at z=0.01 to, say, 7 decimals you need to start with >40 decimals precision.

The solution is of course to avoid this cancellation. A straight-forward way of achieving this is to compute the power series around 0.

You could use sympy to obtain the power series:

>>> z = sympy.Symbol('z')
>>> f = -((z**6-210*z**4+4725*z**2-10395)*sympy.sin(z)+(21*z**5-1260*z**3+10395*z)*sympy.cos(z))/z**7
>>> f.nseries(n=20)
z**6/135135 - z**8/4054050 + z**10/275675400 - z**12/31426995600 + z**14/5279735260800 - z**16/1214339109984000 + z**18/364301732995200000 + O(z**20)

For small z a small number of terms appear to be enough for good accuracy.

>>> ply = f.nseries(n=20).removeO().as_poly()
>>> float(ply.subs(z, 0.1))
7.397541093587708e-12

You can export the coefficients for use with numpy.

>>> monoms = np.array(ply.monoms(), dtype=int).ravel()
>>> coeffs = np.array(ply.coeffs(), dtype=float)
>>> 
>>> (np.linspace(-0.1, 0.1, 21)[:, None]**monoms * coeffs).sum(axis=1)
array([7.39754109e-12, 3.93160564e-12, 1.93945374e-12, 8.70461282e-13,
       3.45213317e-13, 1.15615481e-13, 3.03088138e-14, 5.39444356e-15,
       4.73594159e-16, 7.39998273e-18, 0.00000000e+00, 7.39998273e-18,
       4.73594159e-16, 5.39444356e-15, 3.03088138e-14, 1.15615481e-13,
       3.45213317e-13, 8.70461282e-13, 1.93945374e-12, 3.93160564e-12,
       7.39754109e-12])

Thank you, Paul. I tried to do the plot with all the spherical bessel functions up to 13th order comparing my original method (mpmath and np.vectorize, mp.dps=100), your solution (taking 50 terms into account of the power series, as I need to take into account 0>z>15) and also I found the way to use the implemented scipy functions for my purpose (together with np.vectorize). The fastest of the three options turned out to be the last one , followed by yours, which was just 2.5 times slower, being my initial method by far the slowest (31 times slower). — Edi, Feb 16 '18 at 15:02

Is there a faster way of combining functions containing mpf (mpmath floats) with numpy arrays?

1 Answers1