Nan when the SH order reaches 160

[js1019]

Hi!

I like your code very much. However, I obtained NAN while using inverseSHT when the spherical harmonic order reaches 160. What is the upper bound of the order? Can you also compute higher order transform?

Many thanks,
Jia

[HaHeho]

The cause of this limitation is limited numerical precision. The largest word length for the representation of floating point numbers in Matlab is double precision i.e. 64 bit. This is far more than enough to represent reasonable input and output data for the processing. However, it can be a limiting factor in intermediate steps of the processing. In case of Spherical Harmincs, the getSH() function requires the use of the Matlab factorial() operation.

Spherical-Harmonic-Transform/getSH.m

Line 30 in c3c7419

norm = sqrt( (2*n+1)*factorial(n-m) ./ (4*pi*factorial(n+m)) );

This is what the Matlab documentation mentions about the function:
"For double-precision inputs, the result is exact when n is less than or equal to 21. Larger values of n produce a result that has the correct order of magnitude and is accurate for the first 15 digits. This is because double-precision numbers are only accurate up to 15 digits."
According to the code line above, 21 is reached with order 11 (m+n=11+11=22). Of course, numerical errors increase due to the use of many other operations up to the SH coefficients. However, this is negligible up to a certain point. Where that point is very much depends on your application and the actually employed SH order.

In Matlab, the factorial() function yields Inf above an argument of 170. With SHs, this is reached with order 86. This is where your results have literally zero precision. However, in the realm of spherical microphone array processing (the background of this work by Archontis Politis), the numerical precision is never really an issue. There is simply no microphone arrays that yield high enough spatial resolution without detrimental spatial aliasing errors at much larger magnitude than the errors due to numerical precision.

TL;DR:
I think you should be happy about reaching even close to 160. Now I am curious however. What kind of application requires representations at such a high order?

[js1019]

The cause of this limitation is limited numerical precision. The largest word length for the representation of floating point numbers in Matlab is double precision i.e. 64 bit. This is far more than enough to represent reasonable input and output data for the processing. However, it can be a limiting factor in intermediate steps of the processing. In case of Spherical Harmincs, the getSH() function requires the use of the Matlab factorial() operation.

Spherical-Harmonic-Transform/getSH.m

Line 30 in c3c7419

norm = sqrt( (2*n+1)*factorial(n-m) ./ (4*pi*factorial(n+m)) );

This is what the Matlab documentation mentions about the function:
"For double-precision inputs, the result is exact when n is less than or equal to 21. Larger values of n produce a result that has the correct order of magnitude and is accurate for the first 15 digits. This is because double-precision numbers are only accurate up to 15 digits."
According to the code line above, 21 is reached with order 11 (m+n=11+11=22). Of course, numerical errors increase due to the use of many other operations up to the SH coefficients. However, this is negligible up to a certain point. Where that point is very much depends on your application and the actually employed SH order.

In Matlab, the factorial() function yields Inf above an argument of 170. With SHs, this is reached with order 86. This is where your results have literally zero precision. However, in the realm of spherical microphone array processing (the background of this work by Archontis Politis), the numerical precision is never really an issue. There is simply no microphone arrays that yield high enough spatial resolution without detrimental spatial aliasing errors at much larger magnitude than the errors due to numerical precision.

TL;DR:
I think you should be happy about reaching even close to 160. Now I am curious however. What kind of application requires representations at such a high order?

Thanks for the comments! I was playing with the code to reproduce the topography and crust of Mars and Moon. Please see the link as an example. Though the largest degree of the gravity model of the paper is 150, the highest degree used (in the reference) is 300+.

[HaHeho]

Interesting. I have noticed before that Spherical Harmonics are also utilized in topography. Of course one certainly has to deal with much higher data densities in this context.

So how are typical limitations in numerical precision dealt with in that context? AFAIK it is very cumbersome to get any regular programming language to deal with more than double precision.

[js1019]

Interesting. I have noticed before that Spherical Harmonics are also utilized in topography. Of course one certainly has to deal with much higher data densities in this context.

So how are typical limitations in numerical precision dealt with in that context? AFAIK it is very cumbersome to get any regular programming language to deal with more than double precision.

I guess that one may need to compute the spherical harmonics with associated Legendre polynomials and factorials combined.

[Ragoon35]

usually, the lgamma function is used to avoid large factorial division. Like this:
norm = sqrt( ((2n+1)/(4pi)) * exp(lgamma[m - n]) * exp(-lgamma[m + n]) );