ENH: special: Add new implementations of some of the Mathieu functions by steppi · Pull Request #25128 · scipy/scipy · GitHub updates some of the Mathieu functions to use algorithms worked out by Stuart Brorson (@brorson). It is based on functions from the xsf PR: ENH: Updated and refactored Mathieu function implementation (part 1). by steppi · Pull Request #144 · scipy/xsf · GitHub .
Stuart Brorson’s implementations of the Mathieu functions first computes Fourier coefficients which involves finding eigenvectors of a symmetric tridiagonal recurrence matrix, and then computes the Fourier expansion. To Fourier coefficients depend only on the arguments m and q of the Mathieu functions, and in principle can be reused across different values of x. To enable this reuse, I use a stateful functor for the scalar kernels which caches the previous values of the Fourier coefficients.
For the cache to work correctly, the iteration in the ufunc loops must proceed in such a way that m and q vary most slowly and x varies most quickly. Since it is impossible to control the order of ufunc iteration in such a fine grained way with the current NumPy C API, I’ve added a hack scipy.special._ufunc_tools._with_cache_optimization, which broadcasts all arguments together, and transposes the inputs (and output) so that all of the axes along which variables participating in the cached computation have stridge-length 0 (after broadcasting) are transposed to the end. By forcing iteration in C order, we get the exact iteration behavior we want, and the overhead is essentially negligible.
_with_cache_optimization also supports passthrough of all ufunc kwargs. However, since it replaces the ufunc with a wrapper which is a regular function, we lose ufunc methods and attributes like reduce, nin, nout, etc. There is a similar problem for the delegation set up in _support_alternative_backends, but this is at least hidden behind the SCIPY_ARRAY_API feature flag currently. So that these attributes and methods aren’t lost completely, I’ve added and documented a ufunc attribute to get the raw, non-cache-optimized ufunc.
This is a breaking change, but the benefits are immense. Stuart Brorson’s implementation is significantly more accurate and supports a much wider range of paramters. Without the caching however, the ufunc will recompute the fourier coefficients with an expensive (O(m^3) worst case), lapack call for every value of x for fixed m and q. Adding caching essentially turns an O(m^4) algorithm into an O(m^3) algorithm, and can bring 200x or more speedups in practical workflows.
Is it OK to make this breaking change? One argument for it is that these Mathieu function implementations aren’t currently even production grade, but will be after this PR. I like to imagine anyone relying on them while using ufunc methods would be overjoyed with the improvement in quality, and not mind replacing mathieu_sem.whatever with mathieu_sem.ufunc.whatever if needed. If the consensus is that the breaking change is not justified, one could implement ufunc shims that have the desired methods and attributes, implementing them by passing through to the underlying ufunc.
If interested, or if you have an opinion, please join the discussion at ENH: special: Add new implementations of some of the Mathieu functions by steppi · Pull Request #25128 · scipy/scipy · GitHub.