Wavemoth -- Fast spherical harmonic transforms by butterfly matrix compression

We present Wavemoth, an experimental open source code for computing scalar spherical harmonic transforms (SHTs). Such transforms are ubiquitous in astronomical data analysis. Our code performs substantially better than existing publicly available codes due to improvements on two fronts. First, the computational core is made more efficient by using small amounts of precomputed data, as well as paying attention to CPU instruction pipelining and cache usage. Second, Wavemoth makes use of a fast and numerically stable algorithm based on compressing a set of linear operators in a precomputation step. The resulting SHT scales as O(L^2 (log L)^2) for the resolution range of practical interest, where L denotes the spherical harmonic truncation degree. For low and medium-range resolutions, Wavemoth tends to be twice as fast as libpsht, which is the current state of the art implementation for the HEALPix grid. At the resolution of the Planck experiment, L ~ 4000, Wavemoth is between three and six times faster than libpsht, depending on the computer architecture and the required precision. Due to the experimental nature of the project, only spherical harmonic synthesis is currently supported, although adding support or spherical harmonic analysis should be trivial.

💡 Research Summary

The paper introduces Wavemoth, an experimental open‑source library for scalar spherical harmonic transforms (SHTs) that dramatically outperforms existing public codes. The authors focus on two complementary improvements. First, they optimise the computational core for modern CPUs: a small amount of pre‑computed data is stored, memory accesses are arranged to be cache‑friendly, and instruction‑level pipelining is respected. Second, they adopt the fast, numerically stable algorithm of Tygert (2010), which compresses the set of Legendre‑transform operators (the Λ matrices) using a butterfly matrix compression technique.

The butterfly compression works by recursively applying an Interpolative Decomposition (ID) to sub‑blocks of Λ. In each ID step a low‑rank skeleton of k columns is selected (the “skeleton matrix”) and the remaining columns are expressed as linear combinations via an interpolation matrix. By permuting rows and columns after each level, the matrix is factorised into a product of block‑diagonal matrices S₁,…,Sₚ, permutation matrices P₁,…,Pₚ (the “butterflies”), and a residual block‑diagonal matrix R:

  Λ = R Sₚ Pₚ … S₂ P₂ S₁.

The block size (≈150 × 150 in the authors’ experiments) and the number of levels (determined by column width, 64 columns worked well) are chosen empirically. The compression accuracy is tunable: tighter tolerances yield larger S‑matrices but higher numerical fidelity, while looser tolerances give greater speed‑up with modest loss of precision. In practice the authors achieve near‑double‑precision accuracy even at high compression ratios.

With the compressed representation, the Legendre transform (the dominant O(L N₍ring₎) step) can be performed in O(L² (log L)²) operations, rather than the O(L³) cost of the naïve separated‑sum approach. The FFT part of the SHT (the azimuthal sum) remains unchanged and is already optimal. Consequently, the overall synthesis transform scales as O(L² (log L)²) for the resolution range of interest.

Performance benchmarks are carried out on the HEALPix grid, using libpsht as the reference implementation. For low and medium resolutions (L≈500–1500) Wavemoth is roughly twice as fast as libpsht. At the Planck resolution (L≈4000) the speed‑up ranges from three to six times, depending on the hardware (Intel vs. AMD, cache size) and the requested accuracy. The code is currently limited to spherical harmonic synthesis for real‑valued fields; however, because the synthesis operator is linear, adding analysis (the adjoint transform) would be straightforward.

The paper also surveys other fast SHT approaches—FFT‑based, divide‑and‑conquer, local trigonometric expansions (libftsh), and recent divide‑and‑conquer schemes (SphericalKit)—highlighting their numerical stability issues or lack of publicly available implementations. In contrast, Wavemoth’s algorithm is simple to implement, inherently stable, and delivers substantial gains already at L≈2000.

From a software‑engineering perspective, the authors discuss low‑level optimisation: data layout to minimise cache misses, SIMD vectorisation, and multi‑threading via OpenMP. An appendix provides details of the C implementation and the pre‑computed data format.

In summary, Wavemoth demonstrates that butterfly matrix compression of the Legendre operator, combined with careful CPU‑aware coding, yields a fast, accurate, and openly available SHT library. Its O(L² (log L)²) scaling and demonstrated 3–6× speed‑ups at Planck‑scale resolutions make it a compelling tool for modern cosmological data analysis, high‑resolution sky simulations, and any application requiring large‑scale spherical harmonic synthesis.