Instability in Mittag-Leffler function using NumPy

Question

In trying to reproduce the plot on Wolfram MathWorld, and in trying to help with this SO question, I ran into some numerical instability I don't understand:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

def MLf(z, a):
    """Mittag-Leffler function
    """
    k = np.arange(100).reshape(-1, 1)
    E = z**k / gamma(a*k + 1)
    return np.sum(E, axis=0)

x = np.arange(-50, 10, 0.1)

plt.figure(figsize=(10,5))
for i in range(5):
    plt.plot(x, MLf(x, i), label="alpha = "+str(i))
plt.legend()
plt.ylim(-5, 5); plt.xlim(-55, 15); plt.grid()

You can see the instability best in the orange line where a = 1, starting at about x = -35, but there's a problem with a = 0 (blue line) too. Changing the number of terms to sum (i.e. j) changes the x at which the instability occurs.

What's going on? How can I avoid this?

Answer 1

If a=0, the series definition of MLf that you are using only applies when |z|<1. Indeed, when the base z is greater than 1 in absolute value, the powers z**k keep increasing, and the series diverges. Looking at its 100th, or another, partial sum is pointless, those sums have nothing to do with the function outside of the interval -1 < z < 1. Just use the formula 1/(1-z) for the case a=0.

Case a = 1

The function is exp(z) and technically, it is represented by the power series z**k / k! for all z. But for large negative z this power series experiences catastrophic loss of significance: the individual terms are huge, for example, (-40)**40/factorial(40) is over 1e16, but their sum is tiny (exp(-40) is nearly zero). Since 1e16 approaches the limits of double precision, the output becomes dominated by the noise of truncating/rounding operations.

In general, evaluating polynomials by adding c(k) * z**k is not the best thing to do, both from the efficiency and precision standpoints. Horner's scheme is implemented in NumPy already, and using it simplifies the code:

k = np.arange(100)
return np.polynomial.polynomial.polyval(z, 1/gamma(a*k + 1))

However, this is not going to save the series for exp(z), its numeric issues are beyond NumPy.

You could use mpmath for evaluation, gaining in accuracy (mpmath supports arbitrarily high precision of floating point operations) and losing in speed (no compiled code, no vectorization).

Or you could just return exp(z) from MLf when a=1.

Case 0 < a < 1

The series converges, but again with catastrophic precision loss; and now there is no explicit formula to fall back on. The aforementioned mpmath is one option: set really high precision (mp.dps = 50) and hope it's enough to sum the series. The alternative is to look for another way to compute the function.

Looking around, I found the paper "Computation of the Mittag-Leffler function and its derivative" by Rudolf Gorenflo, Joulia Loutchko & Yuri Luchko; 2002. I took formula (23) from it, and used it for negative z and 0 < a < 1.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.integrate import quad

def MLf(z, a):
    """Mittag-Leffler function
    """
    z = np.atleast_1d(z)
    if a == 0:
        return 1/(1 - z)
    elif a == 1:
        return np.exp(z)
    elif a > 1 or all(z > 0):
        k = np.arange(100)
        return np.polynomial.polynomial.polyval(z, 1/gamma(a*k + 1))

    # a helper for tricky case, from Gorenflo, Loutchko & Luchko
    def _MLf(z, a):
        if z < 0:
            f = lambda x: (np.exp(-x*(-z)**(1/a)) * x**(a-1)*np.sin(np.pi*a)
                          / (x**(2*a) + 2*x**a*np.cos(np.pi*a) + 1))
            return 1/np.pi * quad(f, 0, np.inf)[0]
        elif z == 0:
            return 1
        else:
            return MLf(z, a)
    return np.vectorize(_MLf)(z, a)
    

x = np.arange(-50, 10, 0.1)

plt.figure(figsize=(10,5))
for i in range(1, 5):
    plt.plot(x, MLf(x, i/3), label="alpha = "+str(i/3))
plt.legend()
plt.ylim(-5, 5); plt.xlim(-55, 15); plt.grid()

No numerical issues here.