Efficient evaluation of spherical harmonic expansions

Assume I know that I can express an approximation of a function by $$\sum_{l=0}^{k}\left( \sqrt{A_l} z_{l,0}^1 L_{l,0}(\theta)+\sqrt{2A_l}\sum_{m=1}^l L_{l,m}(\theta)(z_{l,m}^1 \cos(m\phi)+z_{l,m}^2\sin(m\phi))\right)$$ where $$\theta$$ is the inclination and $$\phi$$ azimuth, and $$z_{l,m}^{1,2}$$ are a sequence of indepedent standard normal random numbers, and $$A_l$$ is some pre-decided number so that it decays fast enough. The functions $$L_{l,m}$$ are given by $$L_{l,m}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l,m}(\cos(\theta))$$ where $$P_{l,m}$$ are the associated Lagrange polynomial

of the first type. This I know. I would like to evaluate this sum for

1) Large values of the truncation parameter $$k$$ (like 10^6)

2) many many values of $$\theta$$ and $$\phi$$.

I would like to do this in python, but i suppose I would be ok with writing a compiled program in Fortran for instance, as long as I can acess from python.

My attempt so far is to use scipy.special.lpmv and just blindly evaluate the sum, but that breaks down very fast for not very large values of l,m. There seems to be some sort of trouble with the scipy implementation, see.

https://stackoverflow.com/questions/39249594/numerical-computation-of-the-spherical-harmonics-expansion

import numpy as np
from scipy.special import lpmv sph_harm
from math import factorial
from dolfin import *
from mshr import *
from numpy import sqrt,cos,sin,pi
import numpy as np
no_terms=50
def rand_sequence(noterms,seed):
np.random.seed(seed)
return(np.random.normal(0,1,[noterms,noterms,2]))
def angular_spec(l,alpha):
return((1+l)**(-alpha))
def plm(l,m,theta):
return(sqrt((2*l+1)*factorial(l-m)/4/pi/factorial(l+m))*lpmv(m,l,cos(theta)))    #it seems to break here!!
rands=rand_sequence(no_terms,123)

def GRF(noterms,alpha,rands,x,y,z):
phi=np.arctan2(y,x)
theta=np.arccos(z)
out=0.0
for l in range(noterms):
a=sqrt(angular_spec(l,alpha))*rands[l,0,0]*plm(l,0,theta)
b=sqrt(2*angular_spec(l,alpha))*sum([plm(l,m,theta)(rands[l,m,0]*cos(m*phi)+rands[l,m,1]*sin(m*phi)) for m in range(1,l+1)])

So, any suggestions? Will I have to give up my hope of being able to evaluate for large $$k$$? How can I speed this up. I will have a very large amount of $$\theta$$s and $$\phi$$s to evaluate in...
• Are the $\phi$ values equally spaced? – vibe Oct 29 '19 at 4:27
• The values of $\phi$ are not, sadly. I am essentially evaluating the sums at the vertices of a special tesselation of the sphere. – edo Oct 29 '19 at 10:48
• Your function is a linear combination of normal random variables, so it is a normal random variable itself (mean 0, variance given by a complicated sum of squares of spherical harmonics, probably without a closed form). That might simplify it, but in your code the random variables $z$'s appear to be fixed, being generated once, and not sampled randomly at every evaluation of the function, so maybe not. – Kirill Oct 29 '19 at 16:54