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. – yarchik Apr 24 '18 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. – Jess Riedel Apr 24 '18 at 11:50

• 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[]; ϕ = k[];
x = Sin[ θ] Cos[ϕ];
y = Sin[ θ] Sin[ϕ];
z = Cos[ θ];
{x, y, z}, {k, fMesh // 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[]; ϕ = k[]; w = k[];
w Conjugate[SphericalHarmonicY[10, 9, θ, ϕ] ]
SphericalHarmonicY[10, 9, θ, ϕ], {k, fMesh // 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 = fMesh // N;
Timing[

Norm[

Table[{li, mi} = i; {lj, mj} = j;
Sum[ θ = k[]; ϕ = k[]; w = k[];
w Conjugate[ SphericalHarmonicY[li, mi, θ, ϕ] ]
SphericalHarmonicY[lj, mj, θ, ϕ], {k, m}], {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[]; ϕ = k[]; w = k[];
w Conjugate[SphericalHarmonicY[li, mi, θ, ϕ]],
{k, m}, {i, ord}];

yTab = Table[
{li, mi} = i;
θ = k[]; ϕ = k[]; w = k[];
SphericalHarmonicY[li, mi, θ, ϕ],
{k, m}, {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! – Jess Riedel Apr 24 '18 at 20:40