Analyse too?

[milankl]

CMB.jl/src/sphericalharmonics.jl

Lines 170 to 174 in 2e467aa

	function synthesize_ring(alms::Matrix{C}, θ₀, ϕ₀, nϕ::Integer,
	::Val{pair}) where {C<:Complex, pair}
	R = real(C)
	lmax, mmax = size(alms) .- 1
	nϕr = nϕ ÷ 2 + 1 # real-symmetric FFT's Nyquist length (index)

Hey Justin, many thanks for all your blog articles and AssociatedLegendrePolynomials.jl, great work! We are currently working on SpeedyWeather.jl, a spectral model of the atmosphere. You may have guessed that it also uses a spectral transform, much like you've defined the synthesize_* functions here. We have inherited the numerics from a decade old Fortran model, but want to overhaul the implementation completely. As I have learned loads from your blog articles, especially the one on "Notes on Calculating the Spherical Harmonics" I was wondering whether you have already draft/intent to write a similar one for the analysis, i.e. the inverse transform to the synthesis you've excellently described?

[jmert]

I don't have another article (yet), but I did dust off some notes/code I had lying around and I've pulled together the beginnings of the analysis code. See https://github.com/jmert/CMB.jl/compare/spharms?expand=1.

Note that I haven't written all the tests I intend to yet, so it's very possible that there are some mistakes to be found yet.

[milankl]

Many thanks that's fabulous! I'll have a look at it and will let you know once it's implemented. For SpeedyWeather.jl we need to write a bunch of tests anyway, I'll ping you here when they pass.

[milankl]

Hey Justin, I've used your code as guideline to write a forward and backward spectral transform for SpeedyWeather.jl. This new spectral transform is developed in SpeedyWeather/SpeedyWeather.jl#38 (still looks horrible chaotic at the moment). I got the code all running and it's nice and fast too (60x faster than where we started from, as discussed in SpeedyWeather/SpeedyWeather.jl#32) but we still have some problems with spectral leakage, that means transforming back and forth, doesn't give you (up to rounding errors) exactly the same result.

First thing, I realised that your code assumes equidistant latitudes (θ in your code) which you get from crange/centered_range. Using Gaussian latitudes instead, the relative errors are up to 1% in the spectral leakage. E.g.

julia> alms = zeros(ComplexF64,32,32)
julia> alms[3,1] = 1          # only the l=2,m=0 harmonic
julia> map = gridded(alms)    # uses Gaussian latitudes
julia> alms2 = spectral(map)   
julia> alms2[3,1]             # should be 1+0im
1.0111458144583563 + 0.0im

Using equidistant latitudes this leakage is much reduced

julia> alms2[3,1]
1.0008945580663524 + 0.0im

However, it's still more than one would expect from rounding errors alone. Here's a plot showing absolute error in the coefficients after that back & forth transform

So for l>2 and even and m=0 the error is much larger than it should be. It's similar at different resolutions (here we use lmax=mmax=31 and 96x48 grid points) but most leakage is across degrees not orders and reduces for larger orders (e.g. alms[3,3] = 1 yields a 1e-8 leakage not 1e-3). I was hoping you have some idea where this would come from. I'd appreciate any input in the mean time I'll try to find a solution to that in your blogposts.

[milankl]

Sorry, I've only just understand that _ecp stands for equi-distant cylindrical projection. Makes sense then that it doesn't work as well with Gaussian latitudes! Okay, does that mean in this line of the analysis (aka forward transform, from grid to spectral)

CMB.jl/src/sphericalharmonics.jl

Line 226 in a9907df

a₁, a₂ = (f₁[m+1], f₂[m+1]) .* (sθ * ΔΩ * cis(m * -ϕ₀))

for Gaussian latitudes one multiplies the f1,f2 with the legendre weight of that given latitude instead? Or does that also depend on the normalization? Related, in the synthesis

CMB.jl/src/sphericalharmonics.jl

Lines 171 to 172 in a9907df

	acc₁, acc₂ = (acc₁, acc₂) .* cis(m * ϕ₀)
	i, (acc₁, acc₂) = alias_index(nϕ, m, (acc₁, acc₂))

How would these lines change with Gaussian latitudes? I'm not worried about the aliasing, as we'll run SpeedyWeather.jl always close to the triangular truncation, e.g. lmax=85 and a256x128 grid, hence we're below Nyquist frequencies (if I got this right).

[milankl]

Okay, I think I got it. I changed the sθ * ΔΩ factor to legendre_weight * π/nlat whereby π/nlat is an extra normalization factor on Gaussian latitudes, similar to the ΔΩ = 2π/N_ϕ * π/N_θ that's needed for the equidistant grid. Now the spectral leakage disappeared and all tests pass 🥳

[jmert]

Glad to see you've worked out some of the issues before I've had a chance to respond. I'll still give some comments, though.

1. Yes, ECP means equidistant cylindrical projection, chosen specifically because it's a dead-simple grid.
   1a. I plan on adding a ring-by-ring synthesis/analysis pair of routines, but that's just not implemented yet. That will allow for transformation over any isolatitude (equispaced) ring of pixels because you specify the latitudes explicitly.

2. The choice of sθ * ΔΩ is specifically the naive translation of an integral (\int f(θ,φ) Yℓm^*(θ,φ) sinθ dθ dφ) to a summation, where dΩ = sinθ dθ dφ → sinθ Δθ Δφ → sinθ ΔΩ.

3. The issue of quadrature weights is something I haven't resolved yet, but know it's a problem for ECP and HEALPix grids (which are the two map types I've traditionally used). HEALPix ships their weights, so whenever I get (1a) sorted out that'll help with spectral leakage on HEALPix maps, but that leaves ECP still leaky. The "simplest" brute-force solution is to just run an iteration of analysis+synthesis pairs, e.g.:

function analyze_ecp_iter(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}                                                                                                                                                
    n, m = size(map)                                                                                                                                                                                                                          
    tol = sqrt(eps(maximum(abs, map)))                                                                                                                                                                                                        
    alms = analyze_ecp(map, lmax, mmax)                                                                                                                                                                                                       
    for ii in 2:10  # allow up to 10 iterations before bailing                                                                                                                                                                                
        map′ = synthesize_ecp(alms, n, m)                                                                                                                                                                                                     
        δmap = map - map′                                                                                                                                                                                                                     
        # check for convergence                                                                                                                                                                                                               
        if maximum(abs, δmap) < tol                                                                                                                                                                                                           
            break                                                                                                                                                                                                                             
        end                                                                                                                                                                                                                                   
        alms .+= analyze_ecp(δmap, lmax, mmax)                                                                                                                                                                                                
    end                                                                                                                                                                                                                                       
    return alms                                                                                                                                                                                                                               
end

This is actually an approach well-known in the CMB field — see mentions of iteration in https://healpix.jpl.nasa.gov/html/facilitiesnode7.htm

[milankl]

Many thanks for the clarifications. For our application we are more than happy to just use the regular Gaussian grid. I didn't know you could iteratively increase the accuracy, unfortunately, this is at the moment a bit of an overkill for us. We're more looking towards having a 16-bit transform in the future. As long as the errors are unbiased and only scale with eps, we are (fingers crossed) fine.

On that note, I know from you blog article that you've looked into numerical precision issues for the polynomials, but have you also thought about similar issues within the transform? E.g. one could for example change the order of the summation.

[jmert]

I haven't thought as extensively about numerical issues in the SHT, but I've tried to keep at least an eye on not doing terrible things with the implementation either. For example, using the FFT probably does better than the naive multiplication by exp(imϕ) factors (but I haven't checked). I'm working on a branch (#68) where I'm also going to make use of sincospi and cispi to minimize error from pre-multiplying offsets by a transcendental number (when describing pixelizations on the sphere).

A domain-specific thing to think about is summing from lmax to 0 rather than 0 to lmax — the assumption here is that if your map/spectra are bandwidth limited by something like a beam function, the high-ell terms are smaller in magnitude than the low-ell terms, and if there are sufficiently many small high-ell terms, it's possible for their contributions to get dropped as round-off errors in the sum. In such a case, you can in theory eek out a bit more precision without having to resort to compensated summation by just summing from smallest-to-largest magnitude terms.

FYI, there's now a Part II article on my website describing the analysis steps (which you've obviously already worked out).

[milankl]

think about is summing from lmax to 0 rather than 0 to lmax

Yes, I thought about that. I won't investigate this now, but maybe once we are ready to run everything in Float32, we'll look into this. Reason why I'm asking is that eventually we would like to run the transform in 16-bit arithmetic if that could give us a reasonable speedup. But that obviously depends on many factors...

Many thanks for Part II, I saw it already! Great work and very good to have that as a reference.

[jmert]

If you're considering Float16, I think the most important thing that comes to my mind is to realize that you have to worry about underflow because the exponent ranges get too small for the Legendre polynomial recurrences:

e.g. for a modest 80° latitude:

julia> λlm!(0:100, 0:100, cosd(Float64(10)))
101×101 Matrix{Float64}:
  0.282095   0.0        0.0         0.0          0.0          …   0.0          0.0           0.0          0.0
  0.48118   -0.0599944  0.0         0.0          0.0              0.0          0.0           0.0          0.0
  0.602252  -0.132114   0.0116476   0.0          0.0              0.0          0.0           0.0          0.0
  0.679605  -0.216018   0.0303485  -0.00218464   0.0              0.0          0.0           0.0          0.0
  0.722058  -0.306534   0.0583935  -0.00645435   0.000402371      0.0          0.0           0.0          0.0
  0.733503  -0.399035   0.0961018  -0.0139997    0.00131424   …   0.0          0.0           0.0          0.0
  0.716523  -0.489252   0.143215   -0.0257575    0.00313641       0.0          0.0           0.0          0.0
  0.673486  -0.573237   0.198962   -0.0425756    0.0063136        0.0          0.0           0.0          0.0
  0.606963  -0.647402   0.262107   -0.0651629    0.0113572        0.0          0.0           0.0          0.0
  0.519873  -0.708576   0.331007   -0.094041     0.0188263        0.0          0.0           0.0          0.0
  0.415528  -0.754074   0.403668   -0.129504     0.0293038    …   0.0          0.0           0.0          0.0
  0.297599  -0.781761   0.477823   -0.171585     0.0433682        0.0          0.0           0.0          0.0
  ⋮                                                           ⋱                                           ⋮
 -0.486236  -0.573786   0.559024    0.43218     -0.723674         0.0          0.0           0.0          0.0
 -0.58115   -0.477501   0.641076    0.316862    -0.760539     …   0.0          0.0           0.0          0.0
 -0.658389  -0.366751   0.703937    0.19225     -0.775738         0.0          0.0           0.0          0.0
 -0.715604  -0.24489    0.745721    0.0619967   -0.768805         0.0          0.0           0.0          0.0
 -0.751055  -0.11561    0.765168   -0.0700765   -0.739909         0.0          0.0           0.0          0.0
 -0.763665   0.0171738  0.761688   -0.200088    -0.689847         0.0          0.0           0.0          0.0
 -0.753049   0.149437   0.735381   -0.324212    -0.620025     …   0.0          0.0           0.0          0.0
 -0.719533   0.277171   0.687028   -0.438791    -0.532422         0.0          0.0           0.0          0.0
 -0.664135   0.396506   0.618072   -0.540442    -0.42953         -1.66788e-74  0.0           0.0          0.0
 -0.58854    0.503823   0.530577   -0.626161    -0.314291        -2.30541e-73  2.90362e-75   0.0          0.0
 -0.495047   0.595871   0.427163   -0.693411    -0.190008        -2.25859e-72  4.03383e-74  -5.0548e-76   0.0
 -0.386501   0.669857   0.31093    -0.740196    -0.0602518    …  -1.8109e-71   3.97183e-73  -7.05754e-75  8.79948e-77

julia> λlm!(0:100, 0:100, cosd(Float32(10)))
101×101 Matrix{Float32}:
  0.282095   0.0        0.0         0.0          0.0          …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.48118   -0.0599945  0.0         0.0          0.0              0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.602252  -0.132114   0.0116476   0.0          0.0              0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.679605  -0.216018   0.0303485  -0.00218464   0.0              0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.722058  -0.306534   0.0583936  -0.00645436   0.000402372      0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.733503  -0.399036   0.0961019  -0.0139997    0.00131424   …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.716523  -0.489252   0.143215   -0.0257575    0.00313642       0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.673486  -0.573237   0.198962   -0.0425757    0.00631362       0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.606962  -0.647402   0.262108   -0.065163     0.0113572        0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.519873  -0.708576   0.331007   -0.0940412    0.0188264        0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.415528  -0.754075   0.403668   -0.129504     0.0293039    …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  0.297599  -0.781761   0.477824   -0.171585     0.0433683        0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
  ⋮                                                           ⋱                         ⋮                          ⋮
 -0.486242  -0.573777   0.559027    0.432176    -0.723677         0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.581155  -0.477492   0.641078    0.316857    -0.760541     …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.658393  -0.366741   0.703939    0.192245    -0.775738        -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.715606  -0.244879   0.745721    0.0619909   -0.768804        -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.751056  -0.115598   0.765167   -0.0700823   -0.739907         0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.763664   0.0171857  0.761687   -0.200094    -0.689843         0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 -0.753047   0.149448   0.735379   -0.324218    -0.62002      …   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0
 -0.719529   0.277182   0.687024   -0.438796    -0.532415         0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0
 -0.664129   0.396516   0.618068   -0.540446    -0.429523         0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0
 -0.588533   0.503832   0.530572   -0.626165    -0.314284         0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0
 -0.495039   0.595878   0.427157   -0.693414    -0.19             0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0
 -0.386492   0.669862   0.310923   -0.740198    -0.0602431    …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0

julia> λlm!(0:100, 0:100, cosd(Float16(10)))
101×101 Matrix{Float16}:
 0.282    0.0      0.0       0.0       0.0        0.0        …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.481   -0.05984  0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.6025  -0.1318   0.0116    0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.6797  -0.2157   0.03021  -0.002172  0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.722   -0.3062   0.05814  -0.00642   0.000399   0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.7334  -0.3984   0.0957   -0.01393   0.001304  -7.254e-5   …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.717   -0.4888   0.1426   -0.02562   0.00311   -0.0002575      0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.6743  -0.573    0.1982   -0.04233   0.006264  -0.0006657      0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.6084  -0.648    0.2612   -0.06476   0.01126   -0.001445       0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.522   -0.7095   0.33     -0.0935    0.01866   -0.002789       0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.418   -0.754    0.4026   -0.1289    0.02905   -0.004944   …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.3003  -0.7803   0.4763   -0.1709    0.04297   -0.00821        0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 ⋮                                                ⋮          ⋱                         ⋮                          ⋮
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0        …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0           -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0           -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0        …   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0   0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0            0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0
 0.0      0.0      0.0       0.0       0.0        0.0        …   0.0  0.0   0.0  0.0   0.0  0.0   0.0  0.0  -0.0  0.0

See also jmert/AssociatedLegendrePolynomials.jl#27

[milankl]

Closing this as both synthesis/analysis are available.

[jmert]

Should be closed when #68 is merged.