Bugs in Mathieu functions

[brorson]

This is a known issue. Since the old repo Scipy for the special functions has been closed and this one is now active I wanted to make sure it did not get lost.

The existing Scipy implementation of the Mathieu functions comes from the book "Computation of Special Functions" by Zhang and Jin (their Fortran code has been ported to C by the Scipy maintainers, but the algos remain the same). One step in the Zhang & Jin impl involves finding the roots of a difficult function. The function is defined using recursion, and different roots correspond to different function orders. Since the function is difficult, when looking for the root corresponding to the nth order function, the rootfinder can wander away from the desired root and find an adjacent one which corresponds to the n+1st or n-1st order.

The problem with the Zhang & Jin algo is shown in the below figure, which shows the Mathieu "characteristic values" (eigenvalues) vs. frequency parameter q. You can see jumps from one branch to another one for certain orders and at certain q values.

Compare this with the plot shown on the DLMF (albeit for smaller range)

http://dlmf.nist.gov/28.2.F1

The result is that the actual function can randomly jump from one order to another while the q parameter is smoothly swept over a small range per the below figure.

This is manifestly bad since the code will randomly return the wrong function to the user. The correct behavior is that the Mathieu impl doesn't randomly jump to the wrong branch.

I attach code which reproduces these plots to this bug report below.

What I can offer is an alternate implementation which uses a different method to compute the Mathieu eigenvalues -- one which doesn't jump around since it creates a recursion matrix and then computes its eigenspectrum. I had some undergrads work on this project earlier in the summer and they produced this paper describing the method, as well as some preliminary impls:

https://drive.google.com/file/d/1h6bMnNIUY3zBfP9Mh47vgaA-_Ntc29pZ/view

After their work I wrote up a prototype implementation in Matlab, and included a suite of tests. If there is interest on the part of the Scipy team I can translate the impl to C and submit a patch against the existing C codebase. Before I do that, however, I wanted to file this bug report as a way to kick off the process and gauge any interest on the part of the Scipy team.

Best regards,
Stuart
Northeastern University
Boston, Mass, USA

====================================================

import numpy as np
import scipy as sp
from scipy.special import *
import matplotlib.pyplot as plt

# This fcn creates the plots of ce and the eigenvalues
# evidencing the bug in the Mathieu impls.


N = 200                         # Number of sample pts
v = np.linspace(-180, 180, N)   # Fcn domain -- note this is degrees, not rads. 


#--------------------------------------------------------
# First make plot of eigenvalues vs. q.  This is an expansion
# of the famous plot DLMF Fig. 28.2.1 to larger m and q values.
# The bug is evidenced by the places where the eigenvalues jump
# from one branch of the plot to another.  This happens when the
# rootfinder jumps to the wrong root.
q = np.linspace(0,50,2500)  # Freq parameter

# Eigenvalues for ce
for m in range(9):
    a = mathieu_a(m,q)
    plt.plot(q, a,'b-')

# Eigenvalues for se
for m in range(1,10):
    b = mathieu_b(m,q)
    plt.plot(q, b,'r--')
plt.title("First Matheiu eigenvalues vs. q")
plt.xlabel("q")
plt.ylabel("Eigenvalues")
    
plt.show()

#--------------------------------------------------------
# Now create a figure and 10 vertically stacked subplots.
# For some values of q we get the wrong Mathieu function.
fig, axs = plt.subplots(10, 1, figsize=(8, 12), sharex=True)

# Array of q (frequency parameter) values.
qs = np.array([25.0, 25.1, 25.2, 25.3, 25.4, 25.5, 25.6, 25.7, 25.8, 25.9])
N = len(qs)
m = 6   # Order of Mathieu fcn

for i in range(N):
    q = qs[i]
    ce,_ = mathieu_cem(m,q,v)
    axs[i].plot(v, ce)
    axs[i].grid(True)

axs[-1].set_xlabel('v')
for i in range(N):
    ax = axs[i]
    ax.text(-20,-0.6,'q = %3.1f' % qs[i])        
    
fig.suptitle("Mathieu ce(%d,q,v) for varying q" % m)
plt.show()

[lucascolley]

is this roughly the same issue as scipy/scipy#12267?

[steppi]

Thank you @brorson. That would be much appreciated. How long do you think it will be before you have a patch ready? I'm pretty busy at the moment, but hopefully will get through this busy patch within a few weeks. I should also try to find time to work on adding Arnie Lee Van Buren implementations from here: https://github.com/MathieuandSpheroidalWaveFunctions as reference implementations for the xsf test suite, since currently we are (self) testing against the broken Zhang and Jin implementations.

[brorson]

Thanks for the reply, @lucascolley. Yes, the bug you mention has the same root cause as the bug I report. It also is the source for these other Scipy bug reports:

scipy/scipy#4479
scipy/scipy#14526
https://stackoverflow.com/questions/28205127/mathieu-characteristics-cross-when-plotted

And here is a discussion from last year about the problem where I chime in at the end of the thread:

scipy/scipy#14577

The root cause is that Zhang and Jin used the wrong approach to compute the Fourier coefficients used in the sum. The right thing to do in modern computing is an eigenvalue decomposition. Briscombe, Corless, et al. hint at this in their interesting review paper:

https://arxiv.org/abs/2008.01812

Thank you @steppi. I would guess that I will have something ready to deliver in a month or two. I am also busy with work, and working this project is more of a hobby so the time I can spend on it is constrained. My plan is to finish up the Matlab impl & test suite, then do some testing of the existing Scipy impls using my test suite (just to compare the two), then assuming everything works out go ahead and do the port.

Stuart

[brorson]

By the way, @steppi , I did look at the paper by Van Buren & Boisvert since you suggested it some time ago. Thank for your suggestion -- the paper was very useful and I appreciate the reference. I won't say I grasped every detail in it, but I saw two ideas in the paper: one I liked and one I didn't. But I could be confused and so remain open discussion. The paper I read is this one:

Van Buren & Boisvert

The first idea is presented in section 2. (This is the one I disagree with.) They have an improved method to find roots of the nasty function which gives the ratios of the Fourier coefficients. The nasty function being computed is described in detail in the Briscombe paper above in section 3.4.1. To find the ratio of the Fourier coefficients you do a root-finding procedure on a nasty function which is built using a continued fraction expansion. That function has little, tiny features where the roots are, but is otherwise featureless. That makes it easy for the rootfinder to wander away from the desired root. Fixing this involves making sure the initial guess for the root is very close to the actual root. My own thought is that that doesn't generalize well in the context of libraries for a general-purpose math system like Scipy. It's just too fiddly. Doing an eigenvalue decomposition is very cut-and-dried, so works better for general-purpose applications. But I remain open to the possibility that I am confused about their algo or it doesn't require a lot of fiddling.

Their second idea occurs when computing the radial functions. The radial functions are computed using sums of products of Bessel functions of the same or adjacent order. Here's the expression:

https://dlmf.nist.gov/28.24

equations 2 -- 5. This relation is apparently useful for products of Bessel functions separated by more than one order. The Van Buren & Boisvert idea is to tune the calculation for maximum accuracy by selecting the best separation factor depending upon the inputs to the function. They talk about this in the text below Fig. 5 of the paper. A similar tactic was taken by Bibby & Peterson in their book:

https://link.springer.com/book/10.1007/978-3-031-01717-9

What the REU students have done is nowhere near this level of finesse -- they just used the EVD to get the Fourier coeffs and then used the Bessel product series to compute the radial functions.

Regarding accuracy, my own tests check function values for the angular functions over the entire [-pi,pi] input domain, for q in [-1000, 1000] and n = 0, 1, 2, ... 35. My testing tolerance is around 1e-12, and the angular fcn impls pass all tests except orthogonality for the largest q values where things degrade to errors on order 1e-10.

My suggested next steps are:

1. I need to write up some sort of spec, which lays out the domain over which we have tested the impls, and what accuracy we require over that domain. Unfortunately, most of the time such a spec never exists for any math package, and most people assume the answers they get out of the package are exactly true for all inputs -- a bad assumption.

2. I want to port my tests to Python, then run them against the existing Scipy impls to see how well the existing impl works (when it isn't totally off). It may be that the existing impl works over a larger input domain than ours (or might not), but as far as I know it has never been actually tested for accuracy or domain size. It would make sense for us to at least meet what the existing impl can do (when it isn't totally wrong).

Best regards,
Stuart

[brorson]

Sorry for repeatedly updating this thread. But one final word about the Van Buren & Boisvert strategy for the radial functions -- I believe they propose their method to mitigate the errors which build up when summing a series having alternating plus and minus terms. However, one can do this more simply by just summing the minus and plus terms separately, then adding the two results to get the final sum. I haven't compared their strategy against this simple method (indeed, I only half-way understand their work), but given the choice between a complicated solution and a simple one (which I can understand) I typically go for the simple one first and then do the complicated one next if the simple one doesn't work out.

[steppi]

@brorson, I haven't had time to look into your comment in depth, but please feel free to update the thread as often as you'd like with any insights you might have. I do know that splitting alternating sums into two sums and subtracting at the end is typically not numerically stable, and can be vulnerable to catastrophic cancellation and rounding error. Usually you'd want to do some sort of convergence acceleration technique. I will set up test cases using Van Buren's quad precision implementations as a reference though, so you can have more data to test against.