Package for fast spherical harmonic transform in Mathematica?

Question

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.)

Answer 1

• 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}