pylevin: Efficient numerical integration of integrals containing up to three Bessel functions

The paper introduces pylevin, a Python package that utilizes Levin's method to efficiently and accurately compute highly oscillatory integrals containing up to three Bessel functions, offering performance comparable to specialized single-function tools and significantly outperforming standard quadrature methods for multi-function integrals.

The Gist

Imagine you are trying to listen to a specific, very faint whisper in a room that is absolutely filled with the deafening, rapid-fire sound of a thousand jackhammers. That is essentially what mathematicians and physicists face when they try to calculate certain types of equations involving Bessel functions.

These functions are like the "jackhammer sounds" of the math world. They oscillate (vibrate back and forth) incredibly fast. If you try to measure them using standard tools (like a ruler or a basic calculator), you'll miss the details, get the wrong answer, or spend so much time measuring that you'll never finish.

Here is a simple breakdown of what the paper pylevin is about, using some everyday analogies:

1. The Problem: The "Jackhammer" Noise

In physics, especially when studying things with circular symmetry (like ripples in a pond, sound waves, or the distribution of galaxies in the universe), you often need to add up (integrate) these fast-vibrating Bessel functions.

• The Old Way: Imagine trying to count every single vibration of a hummingbird's wing by looking at it with your naked eye. You might guess, but you'll likely miss thousands of flaps. Standard math tools are like that naked eye; they are too slow and clumsy for these fast vibrations.
• The Existing Solutions: Scientists have built specialized "high-speed cameras" (software packages like FFTLog or hankel) to catch these vibrations. But these cameras are very specific. One camera only works if there is one hummingbird. Another only works if there are two. If you have a flock of three, or if the birds are flying in a weird pattern, those cameras break or can't be used.

2. The Solution: The "Swiss Army Knife" Camera

Enter pylevin. The authors created a new tool that acts like a Swiss Army Knife for these calculations.

• Versatility: Unlike the specialized cameras that only handle one or two birds, pylevin can handle integrals with one, two, or even three Bessel functions at the same time. It doesn't care how complex the pattern is.
• The Secret Sauce (Levin's Method): Instead of trying to count every single vibration (which takes forever), pylevin uses a clever trick called Levin's Method.
  • The Analogy: Imagine you need to know the total distance a car traveled while speeding up and slowing down erratically. Instead of checking the speedometer every millisecond, you look at the car's engine and wheels to understand the pattern of its movement. Once you understand the pattern, you can predict the total distance instantly without checking every single moment.
  • pylevin does this mathematically. It figures out the "pattern" of the vibration and solves the equation based on that, skipping the tedious, slow counting.

3. Why It's a Game Changer

The paper compares pylevin to the old "specialized cameras" and the "naked eye" (standard math tools).

• Vs. Specialized Tools: When the task is simple (just one Bessel function), pylevin is just as fast as the specialized tools. It's not slower, but it's much more flexible.
• Vs. Standard Tools: When the task gets hard (two or three Bessel functions), the difference is night and day.
  • The Analogy: If you ask a standard tool to solve a problem with three Bessel functions, it might take 150 seconds (like walking across the country). pylevin solves the same problem in 0.15 seconds (like taking a bullet train). That is 1,000 times faster.

4. The "Re-use" Superpower

One of the coolest features of pylevin is that it remembers its work.

• The Analogy: Imagine you are baking 1,000 cookies. The dough recipe (the complex math part) is the same for all of them, but the chocolate chips (the specific numbers) change slightly for each batch.
  • A normal baker mixes the dough from scratch for every single cookie.
  • pylevin mixes the dough once, keeps it ready, and just swaps in the new chocolate chips for the next cookie. This makes it incredibly fast when you need to run the same calculation many times with slightly different numbers (which is common in scientific research).

Summary

pylevin is a new, open-source software tool that makes it easy and incredibly fast to solve complex math problems involving vibrating waves (Bessel functions).

• It works for simple cases (1 wave) and complex cases (up to 3 waves).
• It is as fast as the best specialized tools for simple cases.
• It is thousands of times faster than standard tools for complex cases.
• It is built for scientists who need to run these calculations over and over again without waiting hours for the computer to finish.

In short: It turns a math problem that used to take hours into a task that takes a fraction of a second, allowing scientists to focus on discovering new things about the universe rather than waiting for their computers to do the math.

Technical Summary

1. Problem Statement

Integrals involving highly oscillatory Bessel functions are ubiquitous in physical systems with rotational symmetry, particularly in cosmology (e.g., calculating power spectra and angular correlations). These integrals generally take the form:
$$I_{\ell_1 \ell_2 \ell_3}(k_1, k_2, k_3) = \int_a^b dx \, f(x) \prod_{i=1}^N J_{\ell_i}(k_i x)$$
where $N \in \{1, 2, 3\}$, $J_{\ell}$ represents spherical or cylindrical Bessel functions, and $f(x)$ is a non-oscillatory function.

Challenges:

• Oscillations: The rapid oscillations of Bessel functions render standard numerical integration techniques (like adaptive quadrature) inefficient, unreliable, or computationally prohibitive.
• Limitations of Existing Tools: Current specialized software (e.g., hankel, pyfftlog, hankl, pyCCL) is highly optimized but typically restricted to integrals containing at most one Bessel function ($N=1$).
• Flexibility Gap: Methods handling $N=2$ or $N=3$ are either non-existent in public packages or rely on assumptions (like separability of the power spectrum) that limit their applicability in complex scenarios (e.g., modified gravity).

2. Methodology

The paper introduces pylevin, a Python package that implements Levin's method to solve these integrals efficiently.

• Core Algorithm: The method casts the oscillatory integral into the solution of an ordinary differential equation (ODE). It constructs a solution over the integration interval $[a, b]$ using collocation points.
• Adaptive Strategy: The implementation is adaptive. It calculates estimates using $n$ collocation points and $n/2$ points. The interval is recursively bisected in regions where the relative error is largest until convergence is achieved.
• Implementation Details:
  • Language: The core algorithm is written in C++ for performance and wrapped in Python using pybind11.
  • Parallelization: Supports OpenMP for multi-threading.
  • Optimization (Homogeneous Solution Reuse): A key innovation is the ability to reuse the "homogeneous solution" of the underlying linear system. In many applications (e.g., Markov Chain Monte Carlo in cosmology), the structure of the integral remains constant while the non-oscillatory function $f(x)$ changes. pylevin precomputes the homogeneous part, allowing subsequent integrations with new $f(x)$ to be performed an order of magnitude faster.
• Scope: It supports arbitrary orders ($\ell$) and arguments ($k$) for up to three Bessel functions ($N=3$).

3. Key Contributions

1. Generalization: Unlike existing tools limited to $N=1$, pylevin handles integrals with $N=1, 2,$ and $3$ Bessel functions within a single, unified framework.
2. Performance Optimization:
   • Significant speedup via C++ implementation and OpenMP parallelization.
   • Homogeneous solution caching: Dramatically accelerates repeated calculations where only the weight function $f(x)$ changes.
3. Usability: Provides a high-level, user-friendly Python interface that abstracts the complexity of Levin's method.
4. Open Source: The code is available on GitHub and PyPI, released under a CC BY 4.0 license.

4. Results and Benchmarks

The author benchmarked pylevin against specialized codes (hankel, hankl, pyfftlog, pyCCL) and standard adaptive quadrature (scipy.integrate.quad).

• Comparison with Specialized $N=1$ Codes:
  • For single Bessel function integrals (Hankel transforms, FFTLog applications), pylevin is competitive.
  • In tests against hankel (Ogata method) and pyfftlog, pylevin achieved comparable runtimes (within a factor of 2).
  • Against hankl and pyCCL (FFTLog based), pylevin was slightly slower (factor of ~2) but offered greater flexibility (no requirement for separable power spectra or specific grid discretization).
• Comparison with Standard Quadrature ($N=2, 3$):
  • For integrals with two Bessel functions, pylevin was ~3 orders of magnitude faster than scipy.integrate.quad (0.01s vs 45s).
  • For integrals with three Bessel functions, pylevin was ~4 orders of magnitude faster (0.15s vs 150s).
  • Stability: Standard adaptive quadrature failed to converge (catastrophic failure) for large argument ranges, whereas pylevin maintained accuracy and stability.
• Accuracy: In all comparisons where benchmarks existed, pylevin results agreed well with specialized methods, with relative differences typically within numerical noise limits.

5. Significance

• Enabling Complex Cosmological Calculations: pylevin allows for the efficient computation of non-Limber projections and angular power spectra involving multiple Bessel functions without relying on the "separability" assumption required by FFTLog-based methods. This is crucial for testing modified gravity theories where the matter power spectrum may not be separable.
• Bridging the Gap: It fills a critical gap in the numerical analysis toolkit by providing a robust, general-purpose solver for multi-Bessel integrals, which were previously difficult or impossible to compute efficiently.
• Efficiency: By leveraging Levin's method and specific optimizations for repeated integrations, it makes high-precision cosmological parameter estimation (e.g., via MCMC) computationally feasible for complex models.

In summary, pylevin represents a significant advancement in numerical integration for oscillatory functions, offering a versatile, fast, and stable solution that outperforms standard methods by orders of magnitude for multi-Bessel integrals while remaining competitive with specialized single-Bessel solvers.