Announcement: Fast Spherical Harmonics


#1

I’m excited to announce a major release of FastTransforms.jl introducing the first fast and backward stable transform between spherical harmonic expansions and bivariate Fourier series! And it’s in Julia first.

Convert between representations as simply as:

julia> Pkg.add("FastTransforms")

julia> using FastTransforms

julia> F = sphrandn(Float64, 975, 975);

julia> G = sph2fourier(F);

julia> H = fourier2sph(G);

julia> norm(F-H)
1.9353272933547886e-13

The Fortran package SHTOOLS claims to convert spherical harmonic expansions of bandwidths 2,800 to function values on the sphere in around 3 minutes. If we pre-compute a plan, the execution runs roughly 18 times faster:

julia> F = sphrandn(Float64, 2800, 2800);

julia> P = plan_sph2fourier(F; sketch = :none);
Pre-computing...100%|███████████████████████████████████████████| Time: 0:01:11

julia> @time G = P*F;
  9.834737 seconds (6 allocations: 119.608 MiB, 0.23% gc time)

julia> @time H = P\G;
  8.142422 seconds (6 allocations: 119.608 MiB, 0.08% gc time)

julia> norm(F-H)
7.186118715996615e-12

(The two results are essentially related by an FFT, which is negligible in comparison.)

The new algorithms are described in the following pre-print:

R. M. Slevinsky. Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series, arXiv:1705.05448, 2017.