Legendre calculations at ℓ and/or m ≳ 2100 lose accuracy

[jmert]

Take the following setup code:

using PyPlot
import CMB.Legendre: LegendreSphereCoeff

function integrate_simpson(f, xs)
    N = length(xs)
    I = copy(f(xs[1]))
    I .+= f(xs[end])
    for x in xs[2:2:end-1]
        I .+= oftype(x, 2)/3 .* f(x)
    end
    for x in xs[3:2:end-1]
        I .+= oftype(x, 4)/3 .* f(x)
    end
    I .*= step(xs)
    return I
end

# _sinc(x) == sinc(x/π) without the extra division
_sinc(x) = iszero(x) ? one(sin(zero(x))) : sin(x) / x

function pixel_window_function(θmin, θmax, φmin, φmax;
        lmax::Int=720, nstep::Int=100, eltype::Type=Float64)

    nstep += mod(nstep, 2) # Forces nstep to be even
    ell = 0:lmax
    zmin, zmax = cos.(convert.(eltype, (θmax, θmin)))
    zs = range(zmin, stop=zmax, length=nstep)
    Δz, Δφ = eltype.((zmax - zmin, φmax - φmin))

    Λfn(z) = let Λtmp = zeros(eltype, lmax+1, lmax+1),
                 Λtab! = LegendreSphereCoeff{eltype}(lmax)
        Λtab!(Λtmp, z)
    end
    Λ = integrate_simpson(Λfn, zs)
    Λ[:,2:end] .*= sqrt(2) # double to account for m ≤ -1

    Φ(m) = (Δφ * _sinc(m * Δφ/2))^2

    wlm2 = squeeze(sum(Λ.^2 .* Φ.(ell'), dims=2), dims=2)
    norm = 4π / (Δφ * Δz)^2
    Wl = sqrt.(norm ./ (2 .* ell .+ 1) .* wlm2)
    return (ell, Wl, Λ)
end

Computing up to ℓ_max = 2160 for the pixel settings as follow, the answer starts diverging somwhere in the range ℓ ∈ [2120, 2125].

Δδ, Δλ = (0.250,   0.466)
δ₀, λ₀ = (-57.625, 0.000)

θlim = deg2rad.(90.0 .- (δ₀ .+ (Δδ .* (1, -1))))
φlim = deg2rad.(λ₀ .+ (Δλ .* (-1, 1)))
(l, Wl, Λ) = pixel_window_function(θlim..., φlim..., lmax=3*720)

semilogy(l, Wl)

Looking at the Legendre integral term Λ directly, for ℓ = 2110 and ℓ = 2125 (x-axis being across m) we can see what is probably contributing to the divergence as a hump beyond the final "ring-down".

plot(l, Λ[2110, :], label="ℓ = 2110")
plot(l, Λ[2125, :], label="ℓ = 2125")
legend()
xlabel(raw"$m$")
ylabel(raw"$\int\lambda_\ell^m(z)\,dz$")

Repeating the same exercises with eltype = BigFloat do not show the same behavior, so the reasonable conclusion is that there's either round-off or term-by-term cancellation erros which are contributing to the breakdown for double-float precisions.

Other observations:

• Switching from Simpson's rule to trapezoidal rule integration has no impact, so I'm interpreting this as numerical effects from the Legendre polynomials, not overshoot/"ringing" induced by the integration method.

[arturgower]

I imagine this issue would be solved if the package used, say, the associated Legendre polynomial implementation in GSL.jl. Why did you choose to implement this in-house? I am also in need of associated Legendre functions.

[jmert]

I reimplemented because this was initially a proving ground for the method which I then translated into straight C code for use in a Matlab MEX function (that I unfortunately have to use for my research). Compiling in dependencies is really inconvenient in our case (where we're still using Matlab R2009a), so it's better to not have dependences at all.

The other benefit is that once you have everything in Julia, you can do the Really Cool Things™ that Julia empowers — for example, I recently proved to myself that I could get the derivatives of relatively complex functions that themselves depend on the Legendre polynomials by simply leveraging the dual number type from ForwardDiff.jl. That'd be a non-trivial transformation to code up if the Legendre polynomials are coming from a C library like GSL.

I'd love to eventually get back to this and show such a use case, but for the moment my only priority in life is to complete my thesis and graduate ;-)

[arturgower]

Good luck with thesis @jmert, and I totally agree on having as much as possible in Julia.

Maybe one day we both work on getting a stable version of everything spherical harmonics into Julia SpecialFunctions? =]

[jmert]

Note to self: Updating for modern Julia versions and assuming #16 is finished, the pixel window function calculation is instead

using PyPlot
import CMB.Legendre: LegendreSphereCoeff

function integrate_simpson(A::AbstractArray{<:Any,3}, δ)
    N = size(A, 1)
    I = A[1,:,:] .+ A[N,:,:]
    for ii in N-1:-1:2
        @. I = fma(A[ii,:,:], 2(mod(ii,2)+1)/3, I)
    end
    I .*= δ
    return I
end

# _sinc(x) == sinc(x/π) without the extra division
_sinc(x) = iszero(x) ? one(sin(zero(x))) : sin(x) / x

function pixel_window_function(θmin, θmax, φmin, φmax;
        lmax::Int=720, nstep::Int=100, elt::Type=Float64)

    nstep += mod(nstep, 2) # Forces nstep to be even
    ell = 0:lmax
    zmin, zmax = cos.(convert.(elt, (θmax, θmin)))
    zs = range(zmin, stop=zmax, length=nstep)
    Δz, Δφ = (zmax - zmin, φmax - φmin)

    Λtab! = LegendreSphereCoeff{elt}(lmax)
    Λ = zeros(elt, nstep, lmax+1, lmax+1)

    # !!! Need to drop subnormal terms to avoid blow-up at large (l,m)!
    set_zero_subnormals(true)
    Λtab!(Λ, zs)
    set_zero_subnormals(false)

    Λ2 = integrate_simpson(Λ, step(zs))
    Λ2[:,2:end] .*= sqrt(2) # double to account for m ≤ -1

    Φ(m) = (Δφ * _sinc(m * Δφ/2))^2

    wlm2 = dropdims(sum(Λ2.^2 .* Φ.(ell'), dims=2), dims=2)
    norm = 4π / (Δφ * Δz)^2
    Wl = sqrt.(norm ./ (2 .* ell .+ 1) .* wlm2)
    return (ell, Wl, Λ)
end

Δδ, Δλ = (0.250,   0.466)
δ₀, λ₀ = (-57.625, 0.000)

θlim = deg2rad.(90.0 .- (δ₀ .+ (Δδ .* (1, -1))))
φlim = deg2rad.(λ₀ .+ (Δλ .* (-1, 1)))
(l, Wl, Λ) = pixel_window_function(θlim..., φlim..., lmax=3*720)

clf(); semilogy(l, Wl)

Note the calls to set_zero_subnormals() surrounding the Legendre poly calculations — zeroing subnormals stops the blow-up from occurring. On a more fine-grained attempt, explicitly checking for and zeroing subnormals output from just _1term_raise_lm was also sufficient to avoid the accumulated error.