Package for fast spherical harmonic transform in Mathematica?

Mathematica's implementation of the Fast Fourier Transform is, naturally, much faster than computing the discrete transform yourself using Sum.

The analog of the Fourier transform of a function f[theta, phi] on the unit sphere is an expansion in terms of spherical harmonics:

fHarmonics := Table[
NIntegrate[
Conjugate[SphericalHarmonicY[l, m, theta, phi]] f[theta, phi] Sin[theta],
{theta, 0, pi}, {phi, 0, 2 pi}
],
{m, -l, l}, {l, 0, lmax}
]


However, this is very slow, unsurprisingly. Is there a package for Mathematica implementing the widely known algorithms (ref 1, ref 2, ref 3, ref 4) for performing this transformation efficiently?

This was asked without answer in 2013 on Google Groups. Jens Keiner's thesis claims to have implemented this in Mathematica, but I couldn't dig up the code. (I've emailed Keiner and will update this question if he answers.)

• NIntegrate is not so slow by itself, rather the bottleneck is in repeated evaluation of the function for each spherical harmonic. Therefore, a very simple solution is to precompute the function and the weights on a certain mesh once and for all. I typically use the Fibonacci integration (J H Hannay and J F Nye 2004 J. Phys. A: Math. Gen. 37 11591) which is extremely easy to implement and which gives excellent accuracy if the order of SphericalHarmonicY is not too high. Thus, I advocate simplicity (no extra package is needed), although other methods could potentially be much faster. Commented Apr 24, 2018 at 6:38
• @yarchik: Yes, if the function f is expensive to evaluate then pre-computation is useful. (I had already implemented this using FunctionInterpolation on my machine, but didn't want to distract from the main question.) However, even if we ignore the cost of computing f, the integration is slow by itself. Likewise, the FFT is known to be much faster than computing a discrete Fourier transform using Sum even when access to the discrete function values is free; it is the addition and multiplication steps that are slow. Commented Apr 24, 2018 at 11:50
• FYI It seems one of the authors of a nonequispaced fast Fourier transform library is the same Jens Keiner. The lib is written in C, with some M. notebooks in tests folder. Commented Jun 25, 2022 at 9:53

• I am aware that my answer would not be accepted because OP explicitly demanded a FFT-like method.
• I am aware of the fact that the method is not very fast either.
• However, it is so simple that I do this post to help benchmarking more advanced algorithms.
• I follow here the Fibonacci integration algorithm proposed by J. H. Hannay and J. F. Nye, J. Phys. A: Math. Gen. 37 11591 (2004).

Mesh

At first, let us define a function to generate the mesh and corresponding weights. The accuracy depends on the number of points.

fMesh[n_] := Module[{fp, fn, dz, zj, θ, ϕ, w},
fp = Fibonacci[n - 1];
fn = Fibonacci[n];
dz = 2/fn;
Flatten[Table[
zj = -1 + j dz;
θ = ArcCos[zj + Sin[π zj]/π];
ϕ = π j fp/fn;
w = π dz (1 + Cos[π zj]);
{{θ, ϕ, w}, {θ, -ϕ, w}}, {j, 0, fn}], 1]
]


Plot a representative mesh (19th order)

v = Table[θ = k[[1]]; ϕ = k[[2]];
x = Sin[ θ] Cos[ϕ];
y = Sin[ θ] Sin[ϕ];
z = Cos[ θ];
{x, y, z}, {k, fMesh[19] // N}];
ListPointPlot3D[v, BoxRatios -> 1, ColorFunction -> "Rainbow",
Boxed -> False, Axes -> None]


Simplest application---normalization of $$Y_{10}^9(\theta,\phi)$$

Verify the normalization of a particular spherical harmonics

Sum[ θ = k[[1]]; ϕ = k[[2]]; w = k[[3]];
w Conjugate[SphericalHarmonicY[10, 9, θ, ϕ] ]
SphericalHarmonicY[10, 9, θ, ϕ], {k, fMesh[15] // N}]


Brute-force evaluation of many integrals

Use a brute-force approach to compute the integrals of two spherical harmonics. The output should be an identity matrix. The norm of its deviation from the Identity matrix will be our error measure.

lmx = 4;
ord = Flatten[Table[{l, m}, {l, 0, lmx}, {m, -l, l}], 1];
m[15] = fMesh[15] // N;
Timing[

Norm[

Table[{li, mi} = i; {lj, mj} = j;
Sum[ θ = k[[1]]; ϕ = k[[2]]; w = k[[3]];
w Conjugate[ SphericalHarmonicY[li, mi, θ, ϕ] ]
SphericalHarmonicY[lj, mj, θ, ϕ], {k, m[15]}], {i, ord}, {j, ord}]

-IdentityMatrix[Length[ord]]
]
]


lmx-maximal angular momentum, ord-list of {l,m} values. The computation takes around 4 seconds and yields quite precise results.

{4.64707, 1.02165*10^-9}

Matrix multiplication for speed

Now we do the same, however, we use the matrix multiplication in order to compute all the integrals at once.

Timing[

fTab = Table[
{li, mi} = i;
θ = k[[1]]; ϕ = k[[2]]; w = k[[3]];
w Conjugate[SphericalHarmonicY[li, mi, θ, ϕ]],
{k, m[15]}, {i, ord}];

yTab = Table[
{li, mi} = i;
θ = k[[1]]; ϕ = k[[2]]; w = k[[3]];
SphericalHarmonicY[li, mi, θ, ϕ],
{k, m[15]}, {i, ord}];

Norm[Transpose[fTab].yTab - IdentityMatrix[Length[ord]]]
]


The result is identical, however, the computation is >10 times faster

{0.257956, 1.02165*10^-9}

• Useful, thanks! Commented Apr 24, 2018 at 20:40
• It works well. Thanks. Commented Jun 1, 2020 at 14:56
• @tanghe2014 You are welcome, happy it was helpful. Commented Jun 1, 2020 at 15:01
• @yarchik The mesh duplication implemented here differs from the one described by Hannay & Nye. On the line {{θ, ϕ, w}, {θ, -ϕ, w}}, the duplication by Hannay & Nye puts the second mesh at {{θ, ϕ, w}, {θ, ϕ + π, w}}, i.e., rotated by $\pi$ from the first one about the $z$ axis, instead of reflected on the $xz$ plane. The reflected version adds some anisotropy to the mesh; this isn't terrible (as shown by your benchmarks), but it does open the question on whether the altered duplication degrades the performance for integrands with higher anisotropy than the $Y_{lm}$s. Commented Jun 21, 2021 at 15:57
• @EmilioPisanty Interesting observation. I know that you posted a question on this topic. Maybe you can answer your own question and demonstrate the superiority of the Hannay and Nye choice. Commented Jun 21, 2021 at 17:20