Stable Hermite transforms via the Golub-Welsch algorithm

We introduce an efficient and stable algorithm for transforms associated with expansions in Hermite functions interpolated at Hermite polynomial roots. The Hermite transform matrix can be factorised into a diagonal component and an orthogonal matrix, leading to a form which allows both the forward and inverse Hermite transforms to be computed stably. Our novel algorithm computes this factorisation based on the eigendecomposition of the Jacobi matrix associated with Hermite functions. Through numerical experiments, we demonstrate the stability and efficiency gains of this novel method over prior work. Numerical experiments show that the new approach matches or improves on the accuracy of existing stabilized methods, is substantially faster in practice, and enables reliable use of large Hermite expansions in downstream PDE computations. We also provide an open-source implementation, together with reference implementations of previous methods, to facilitate adoption by the community.