Efficiently generating 2-D Gaussian random fields on the sphere

Question

Context

A couple of years ago I posted this question for an efficient code to generate an n-D Gaussian random field (sometimes called processes in other fields of research), which has applications in cosmology.

For instance, the function GaussianRandomField would work as follows in 2D

 u = GaussianRandomField[] //Chop// GaussianFilter[#, 1] &;
 Image[u] // ImageAdjust

Question

I am interested in an efficient way to proceed on the sphere.

Motivation

I would eventually like to make maps like this map of the Cosmic Microwave Background seen by the Planck satellite.

using this code.

Possible difficulties

It most likely involves fast Spherical Harmonic transforms of relatively high order? One possibility might be to link to the Healpix library, but hopefully it would be an overkill.

Feeble attempt

Extracting the transformation code from the above mentioned page

invmollweide[{x_, y_}] := 
 With[{theta = ArcSin[y]}, {Pi x/(2 Cos[theta]), 
   ArcSin[(2 theta + Sin[2 theta])/Pi]}]

I can generate a map:

     Clear[alms]; lmax = 48;
Do[alms[l,  m] = (Random[NormalDistribution[]] + I Random[NormalDistribution[]])
  /Sqrt[(l + 2) (l + 1)]; 
  alms[l, -m] = (-1)^m  Conjugate@alms[l, m];, {l, 0, lmax}, {m, 0, l}];

Do[alms[l, 0] = (Random[NormalDistribution[]])/Sqrt[(l + 2) (l + 1)]; , {l, 0, lmax}];

Clear[field];
field[θ_, ϕ_] = Sum[alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l}];

fieldN = Compile[{θ, ϕ}, field[θ, ϕ] // Evaluate];

dat = ParallelTable[fieldN[θ, ϕ], {θ, 0, Pi, Pi/128.}, {ϕ, 0.,2 Pi, 2 Pi/256.}] // N;

and then plot it:

im = Re[dat] // Image // ImageAdjust// Colorize[#,
 ColorFunction -> "LightTemperatureMap"] &;
 mol=ImageTransformation[im, invmollweide, 
 DataRange -> {{-Pi, Pi}, {-Pi/2, Pi/2}}, 
 PlotRange -> {{-2, 2}, {-1, 1}},Padding-> White]

But

1. The production of the map is rather slow…

2. The code breaks down e.g. at lmax=64

probably because of the accuracy of the spherical harmonics (whereas in astronomy people routinely use lmax = 4096 or more).

It's a pity because Mathematica allows for some cool visualization.

SphericalPlot3D[1, {u, 0, Pi}, {v, 0, 2 Pi}, PlotPoints -> 50, 
 MaxRecursion -> 0, Mesh -> True, 
 TextureCoordinateFunction -> ({#5, 1 - #4} &), 
 PlotStyle -> Directive[Texture[im], Specularity[White, 50]], 
 Lighting -> "Neutral", Boxed -> False, Axes -> False]

See also this

Possible Extension: random vector field on sphere

One can also define vector fields on the sphere as follows

fc[phi_] := 
  Block[{theta}, 
   If[Abs[phi] == Pi/2, phi, theta /. 
     FindRoot[2 theta + Sin[2 theta] == Pi Sin[phi], {theta, phi}]]];
cart[{lambda_, phi_}] := 
 With[{theta = fc[phi]}, {2/Pi*lambda Cos[theta], Sin[theta]}]

Then, if we define VectorSphericalHarmonicV

Clear[ϵ];(*Polarization vector*)
ϵ[λ_] = Switch[λ, -1, {1, -I, 0}/Sqrt[2], 0, {0, 0, 1}, 
  1, {1, I, 0}/Sqrt[2]];

Clear[VectorSphericalHarmonicV];
VectorSphericalHarmonicV[ℓ_, J_, M_, θ_, ϕ_] /; 
  J >= 0 && ℓ >= 0 && Abs[J - ℓ] <= 1 &&  Abs[M] <= J := 
 Sum[If[Abs[M - λ] <= ℓ, 
    ClebschGordan[{ℓ, M - λ}, {1, λ}, {J, M}], 0]*
   SphericalHarmonicY[ℓ,   M - λ, θ, ϕ]*ϵ[λ], {λ, -1, 1}]

pp = Sum[(Random[NormalDistribution[]] + I Random[NormalDistribution[]]) 
Rest@ VectorSphericalHarmonicV[i, i - 1, i - 1, θ, ϕ], {i,2, 8}];

We can write

   gr2 = StreamPlot[pp//Re, {ϕ, -Pi, Pi}, {θ, -Pi/2, Pi/2}, 
   AspectRatio -> 1/2, Frame -> False, 
   StreamColorFunction -> "ThermometerColors", StreamPoints -> 250];
   gr2 = gr2 /. Arrow[pts_] :> Arrow[(cart /@ pts)] /. 
    Point[pts_] :> Point[cart[pts]] // 
   Show[#, PlotRange -> {{-2, 2}, {-1, 1}}] &;

and

Graphics[{Inset[mol, {-2, -1}, {0, 0}, {4, 2}], First[gr2]}, 
 PlotRange -> {{-2, 2}, {-1, 1}}]

I guess once we have a fast spherical harmonic code on the sphere it will be trivial to generalize it to vector fields.

Update

This package might be of relevance. I quote

Wolfram Mathematica is a powerful and convenient software package used by many cosmologists everywhere. However, since it is not always the most efficient for low-level computing, most popular algorithms for numerical computations in cosmology have been written in C or Fortran, since these languages are typically much better suited for the task at hand. This package bundles the functionality of some of these algorithms in a Mathematica package, which makes them easier to use and avoids the need to learn C or Fortran.

Answer 1

A part from efficiency, I noticed that with the definition below, lmax >= 64 works:

Clear[field];
field[θ_, ϕ_] :=
    Chop@
       Total[
          Table[
            alms[l, m] SphericalHarmonicY[l, m, θ, ϕ],
            {l, 0, lmax}, {m, -l, l}
          ]
      , 2
    ];


nn = 4.;
dat = ParallelTable[
    field[θ, ϕ],
         {θ, 0, Pi, Pi/nn},
         {ϕ, 0.,  2 Pi, 2 Pi/nn/2}
    ];

Round[Re[dat], 0.01]

{ {0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51},

{0.02, -1.54, -1.48, -0.47, -2.29, 0.61, 2.05, 1.82, 0.02},

{0.24, -1.61, 0.44, 0.32, -0.55, 0.65, -1.12, -0.08, 0.24},

{-0.34, 0.63, -0.54, 2.08, -0.72, 1.09, -2.09, -1.3, -0.34},

{0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} }

which looks better and produces a reasonable plot.

So, I tested also:

Clear[field];
field[θ_, ϕ_] :=
   Sum[
      alms[l, m] SphericalHarmonicY[l, m, θ, ϕ],
      {l, 0, lmax}, {m, -l, l}
   ];

and this too seems to me to work.

Consequently, it appears to me that Compile over Sum is doing some inappropriate manipulation.

EDIT

I could perhaps track down the problem a little better: it seems to be the SphericalHarmonicY at angles θ below Pi/4 and above Pi - Pi/4, which finally provide extremely high values in dat and hides the details of the other values (at least for lmax = 64).

For example the table:

dat = ParallelTable[
    fieldN[θ, ϕ],
        {θ, Pi/4., Pi - Pi/4., Pi/128.},
        {ϕ, 0, 2 Pi, 2 Pi/256.}
    ];

seems to me to produce a reasonable chart.

Answer 2

GaussianRandomField is a special case. More generally, what is required is fast code for the (inverse) Spherical Harmonic Transform (SHT), which will work for any coefficients $a_{l,m}$. SHTns is a high performance library for Spherical Harmonic Transform written in C and so should be straightforward to link in using MathLink. It would be very nice if this code was built-in to a future version of Mathematica.

To see why fast code is required, for $l_{max} = 2^{12}$ direct computation requires evaluation of around $2^{24}$ spherical harmonics on your grid of $\theta$ and $\phi$ values (another factor of $2^{15}$). There are recursive tricks, of course, but fast methods rely on the FFT.

Regarding your example code, I do not see why you have (-1)^m in your definition for alms[l, -m]. To obtain an explicitly real function you should just have alms[l, -m] = Conjugate@alms[l, m].

Further to your visualisation on the sphere, you can also select your view-centre, e.g. Galactic Coordinates.