Stable Hermite transforms via the Golub-Welsch algorithm

1. Problem Statement

The paper addresses a critical numerical instability encountered in Hermite spectral methods used for solving Partial Differential Equations (PDEs), such as the semi-classical Schrödinger equation and the Gross–Pitaevskii equation.

• The Core Issue: Standard algorithms for computing the Hermite transform matrix $T$ (which maps expansion coefficients to function values at Gauss–Hermite nodes) and its inverse $T^{-1}$ suffer from catastrophic underflow and overflow in double-precision arithmetic when the expansion degree $N$ exceeds approximately 766.
• The Cause: The Hermite functions $\psi_n(x)$ involve a Gaussian term $e^{-x^2/2}$. For large $N$, the nodes $x_k$ extend to large magnitudes where $e^{-x^2/2}$ underflows to zero, causing the loss of information in the first columns of the matrix. Conversely, intermediate calculations can overflow.
• Impact: This instability limits the applicability of Hermite spectral methods in problems requiring high resolution (large $N$), such as resolving the mass ratio between electrons and atomic nuclei in quantum mechanics.
• Existing Limitations:
  • Direct Recurrence: Unstable for $N \ge 766$.
  • Fast Transforms ($O(N \log N)$): Unstable for $N \ge 64$.
  • Bunck's Stabilized Recurrence: Stable but computationally slower than the direct method and less widely known.

2. Methodology

The authors propose a novel algorithm based on the Golub–Welsch algorithm, which computes Gaussian quadrature rules via the eigendecomposition of a Jacobi matrix.

Key Mathematical Insight

The paper establishes a factorization of the Hermite transform matrix $T$ into a diagonal matrix $D$ and an orthogonal matrix $Q$:
$$T = DQ^T$$
$$T^{-1} = QD^{-1}$$
Where:

• $Q$ is an orthogonal matrix derived from the eigenvectors of the Hermite Jacobi matrix.
• $D$ is a diagonal matrix containing scaling factors related to the Hermite functions at the nodes.

The Algorithm (Golub–Welsch Approach)

1. Jacobi Matrix Construction: Construct the tridiagonal Jacobi matrix $J$ associated with the Hermite polynomials using the recurrence coefficients.
2. Eigendecomposition: Compute the eigendecomposition $J = Q \Lambda Q^T$.
   • The eigenvalues $\Lambda$ provide the Gauss–Hermite quadrature nodes ($x_k$).
   • The eigenvectors form the orthogonal matrix $Q$.
3. Diagonal Scaling ($D$): Compute the diagonal entries $D_{jj} = \sqrt{N}|\psi_{N-1}(x_j)|$.
   • To avoid underflow/overflow when computing $\psi_{N-1}(x)$, the authors use a hybrid approach:
     • Small $N$ ($N < 200$): Use Clenshaw's algorithm (recurrence).
     • Large $N$ ($N \ge 200$): Use an asymptotic expansion involving parabolic cylinder functions and Airy functions.
4. Reconstruction: The final transform matrices are reconstructed using the computed $Q$ and $D$.

Comparison with Bunck's Method

The paper also provides a modern exposition of Bunck's stabilized recurrence (Algorithm 1), which scales the vectors during computation to keep values within a safe range ($O(1)$) and rescales them at the end. While stable, the authors note it is computationally more expensive than their new approach.

3. Key Contributions

1. Novel Stable Algorithm: Introduction of a stable $O(N^2)$ algorithm for computing Hermite transforms based on the Golub–Welsch framework, avoiding the underflow/overflow issues of direct recurrence.
2. Theoretical Unification: Demonstrating that the Hermite transform matrix factorization is intrinsically linked to the eigendecomposition of the Hermite Jacobi matrix, providing a conceptually simple route to stable computation.
3. Hybrid Evaluation Strategy: Developing a robust method to compute the necessary diagonal scaling factors ($D$) by switching between Clenshaw's algorithm and asymptotic expansions depending on $N$.
4. Open-Source Implementation: Providing a public MATLAB implementation and reference implementations of previous methods to facilitate community adoption.
5. Free Quadrature Data: The method yields Gauss–Hermite nodes and weights "for free" as by-products of the eigendecomposition, whereas other methods often require pre-computation of these values.

4. Results

The authors validate their method through extensive numerical experiments comparing three approaches: Direct, Bunck, and Golub–Welsch.

• Accuracy:
  • The Direct method fails (error $\approx 1$) for $N \ge 766$.
  • Both Bunck and Golub–Welsch maintain high accuracy (machine precision, $\approx 10^{-15}$) even for $N$ up to $1000$.
  • The transition to the asymptotic expansion at $N=200$ in the Golub–Welsch method does not degrade accuracy.
• Efficiency (Speed):
  • The Golub–Welsch method is significantly faster than Bunck's method.
  • It is comparable in speed (and marginally faster in some cases) to the unstable Direct method, offering the best of both worlds: stability and speed.
• Downstream Application (Gross–Pitaevskii Equation):
  • The authors applied the transforms to a splitting method for the Gross–Pitaevskii equation (modeling Bose–Einstein condensation).
  • Simulations using the Direct method crashed or became unstable around $N \approx 766$.
  • Simulations using the Golub–Welsch method remained stable and accurate for $N$ up to 1000, demonstrating its utility for practical, high-resolution PDE simulations.

5. Significance

This work resolves a long-standing barrier in the application of Hermite spectral methods. By enabling stable computations for large $N$ ($>766$), it allows researchers to:

• Solve semi-classical quantum mechanics problems that require extremely high resolution.
• Simulate dispersive PDEs on unbounded domains with high fidelity.
• Utilize Hermite expansions in scenarios where previous "fast" or "direct" methods were numerically unreliable.

The paper effectively bridges the gap between classical Gaussian quadrature theory (Golub–Welsch) and modern spectral method implementation, providing a robust, efficient, and accessible tool for the scientific computing community.