More Notes on Calculating the Spherical Harmonics

Analysis of maps to harmonic coefficients

Notes on Spherical Harmonics Series: Parts 1, 2, 3, 4

This article is a long-overdue follow up to Notes on Calculating the Spherical Harmonics which considers the analysis of real-space map on the sphere to its spherical harmonic coefficients.


Introduction

In Notes on Calculating the Spherical Harmonics, I shared notes on the synthesis of real-space images given a set of spherical harmonic coefficients, but I didn’t discuss the reverse operation of analysis of images back into harmonic coefficients. This article completes the sequence and discusses the analysis tranform.

As a quick reminder, the spherical harmonics \(Y_{\ell m}(\theta,\phi)\) can be written as \begin{align} Y_{\ell m}(\theta, \phi) &\equiv \sqrt{\frac{2\ell+1}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}} P_\ell^m(\cos\theta) \, e^{im\phi} = \lambda_\ell^m(\cos\theta) \, e^{im\phi} \end{align} where \(P_\ell^m(\cos\theta)\) are the Associated Legendre Polynomials, \(\lambda_\ell^m(\cos\theta)\) are defined as the Legendre polynomials with the normalization factor baked in (see my Legendre polynomials series), and \(e^{im\phi}\) are the complex exponentials.

As a complete and orthonormal basis for functions on the sphere, an arbitrary function on the sphere \(f(\theta,\phi)\) can be harmonically decomposed into an ensemble of coefficients \(a_{\ell m}\): \begin{align} a_{\ell m} &= \int_\mathbb{S} f(\theta, \phi) \overline{Y_{\ell m}}(\theta,\phi) \,\sin\theta\,d\theta\,d\phi \end{align} where \(\ell\) and \(m\) are positive integers, \(|m| \le \ell\), and the overline on \(\overline{Y_{\ell m}}\) denotes complex conjugation. The function is resynthesized from the \(a_{\ell m}\)s as the linear combination of the spherical harmonics: \begin{align} f(\theta,\phi) &= \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m}(\theta,\phi) \label{eqn:complex_synthesis} \end{align} Part I described the synthesis step, concentrating on enhancements that go beyond the trivial translation of Eqn. \(\ref{eqn:complex_synthesis}\) to code. In particular:

  1. Symmetries of the sphere permit vastly reducing the computational cost of performing the synthesis transform. Constructing a pixelization which (a) has isolatitude rings of pixels and (b) has north-south mirror symmetry over the sphere enables optimizations that allow computing the (expensive) associated Legendre polynomials for only half the number of rings over the sphere.
  2. Equispaced rings allow recasting the complex exponential summation in \(m\) to be performed via a fast Fourier transform (FFT). Notably, the previous article describes in detail how to properly account for aliasing of modes for \(m > N_\phi\) (where \(N_\phi\) is the number of pixels in the ring).

Implementing these two classes of optimizations resulted in a highly-performant spherical harmonic synthesis algorithm.

We now consider the reverse operation: spherical harmonic analysis.

Spherical Harmonic Analysis

For a discrete image over the point set \((\theta_i, \phi_i)\) (where \(i = \{1, 2, \ldots, N\}\) for \(N\) pixels), the integral must necessarily be discretized: \begin{align} a_{\ell m} = \int_\mathbb{S} f(\theta, \phi) \overline{Y_{\ell m}}(\theta,\phi) \,\sin\theta \,d\theta \,d\phi \qquad \leftrightarrow \qquad \hat a_{\ell m} = \sum_{i=1}^{N} f(\theta_i, \phi_i) \overline{Y_{\ell m}}(\theta_i, \phi_i) \,\sin\theta_i \,\Delta\theta_i \Delta\phi_i \label{eqn:analysis_integral_sum} \end{align} where \(\Delta\theta_i \Delta\phi_i\) is the coordinate area of the \(i\)th pixel. The naive implementation is remarkably similar to the synthesize_reference() function given in Part I with only minor changes necessary. Here we continue to assume that the input map is an equidistant cylindrical projection (ECP) over the sphere.

using AssociatedLegendrePolynomials: λlm!

# Range of pixel centers that cover the interval between `a` and `b` in `n` steps.
function crange(a, b, n)
    r = range(a, b, length = n+1)
    s = step(r)
    return range(a + s/2, b - s/2, length = n)
end

function analyze_reference(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}
    ,  = size(map)
    θ = crange(0.0, 1.0π, )
    ϕ = crange(0.0, 2.0π, )
    ΔΩ = 2R(π)^2 / ( * )

    C = complex(R)
    Λ = Matrix{R}(undef, lmax + 1, mmax + 1)
    Φ = Vector{C}(undef, mmax + 1)
    alms = zeros(C, lmax + 1, mmax + 1)

    for y in 1:, x in 1: 
        ,  = sincos(θ[y])
        λlm!(Λ, lmax, mmax, )       # λ_ℓ^m(cos θ) factors
        Φ .= cis.((0:mmax) .* -ϕ[x])  # e^{-imϕ} factors

        # Σ_{ℓ = 0}^{ℓmax}
        for  in 0:lmax
            # Σ_{m = 0}^{mmax}
            alms[+1,1] += map[y,x] * Λ[+1,1] *  * ΔΩ
            for m in 1:min(, mmax)
                alms[+1,m+1] += map[y,x] * Λ[+1,m+1] * Φ[m+1] *  * ΔΩ
            end
        end
    end
    return alms
end

Like with the reference synthesis implementation, this naive translation of math to code is computationally slow, but it is again a rather direct (and therefore easy to understand) transformation.

To demonstrate the ability to analyze a field, we’ll take the map which was synthesized in Part I and analyze it back to its spherical harmonic coefficients.

A particular realization of a Guassian random field on the sphere, where the modes were drawn from a Guassian distribution and scaled in \(\ell\) such that the variance is proportional to \(\ell^{-2}\).
alms = @time analyze_reference(map, 100);
 20.217139 seconds (2.50 M allocations: 38.382 MiB)

To quantitatively assess whether we’ve achieved a faithful analysis of the given map, we must simply compare the output to the known input. Because the coefficients are complex numbers in (half of) a large square matrix, visually assessing element-by-element is less convenient than summarizing to a one-dimensional, real line, and the obvious choice is to compare the angular power spectra since the input harmonic coefficients were drawn from a distribution consistent with a known spectrum. More specifically, the angular power spectrum \(C_\ell\) for a single map is given by \begin{align} C_\ell = \frac{1}{2\ell + 1} \sum_{m=-\ell}^\ell a_{\ell m} \overline{a_{\ell m}} = \left\langle |a_{\ell m}|^2 \right\rangle \end{align} where the angle brackets \(\langle \cdot \rangle\) denote the sample expectation value. Note that the angular power spectrum is simply a measure of the variance of modes at each \(\ell\) over all available \(m\).

Figure 2.2 below shows the angular power spectra of the input \(a_{\ell m}\) set and those reconstructed by analyzing the synthesized map shown in Figure 2.1.

A comparison of the angular power spectra \(C_\ell\) for the original (“true”) harmonic coefficients (blue line) and for the analyzed map (orange dashed line), both compared to the reference \(C_\ell = \ell^{-2}\) spectrum from which the random gaussian modes were drawn. Because the synthesis and analysis spectra agree very well in the upper panel, the lower panel shows the fractional difference between the two.

We have apparently achieved sub-percent precision in analysis of an image on the sphere into its spherical harmonic components. Before discussing further whether this is as good as can be achieved, though, let us first discuss how to reuse the symmetries discussed at length in Part I to greatly speed up the analysis operation.

Fast Spherical Harmonic Analysis

To achieve much faster analysis, we proceed much like in Part I’s discussion of the fast spherical harmonic synthesis algorithm. Starting from Eqn. \(\ref{eqn:analysis_integral_sum}\), we moderately rewrite the sums to explicitly sum over the independent colatitude \(\theta\) and azimuth \(\phi\) coordinates. (Note that I’ve dropped explicit subscript indices for brevity, but it should not be forgotten that the map is discretized.) \begin{align*} a_{\ell m} &= \sum_{\theta} \sum_{\phi} f(\theta, \phi) \overline{Y_{\ell m}}(\theta, \phi) \,\sin\theta \,\Delta\theta \Delta\phi \\ {} &= \sum_{\theta} \left[ \sum_{\phi} f(\theta, \phi) e^{-im\phi} \right] \lambda_\ell^m(\cos\theta) \,\sin\theta \,\Delta\theta \Delta\phi \end{align*} Again, this form is highly suggestive that a discrete Fourier transform can be employed to transform this double loop with \(\bigO(N^2)\) complexity into a calculation which instead scales as \(\bigO(N\log N)\). We continue to demand that there be isolatitude rings, and that each isolatitude ring of pixels must be uniformly spaced in azimuth so that the sequence of pixel azimuths are \begin{align*} \phi_k = \phi_0 + \frac{2\pi k}{N_\phi} \end{align*} where \(\phi_0\) is the offset of the first pixel in the ring away from \(\phi = 0\), \(N_\phi\) is the number of pixels in the ring, and \(k\) is an integer in \([0, N_\phi - 1]\). Then the expression can be rewritten as \begin{align} a_{\ell m} &= \sum_{\theta} h_m \lambda_\ell^m(\cos\theta) \,\sin\theta \,\Delta\theta &\text{where } h_m \equiv e^{-im\phi_0} \Delta\phi \sum_{\phi} f(\theta, \phi) e^{-2\pi imk/N_\phi} \label{eqn:fast_analysis} \end{align} where we’ve again separated and rewritten the terms such that the factor \(h_m\) is the discrete Fourier transform of the map ring over azimuth, which can be computed via the fast Fourier transform. For the ECP grid being assumed here, the \(\Delta\theta\) and \(\Delta\phi\) are constant factors over the entire sphere and equal to \(\pi/N_\theta\) and \(2\pi/N_\phi\), respectively.

Having isolatitude rings means that as with the synthesis algorithm, it is much more efficient to calculate the associated Legendre polynomials only once per ring rather than for every pixel as the naive implementation above does. Furthermore, we can immediately invoke the discussion from Part I to analyze north/south ring pairs which mirror each other over the equator — the Legendre polynomial terms have the symmetry that \(\lambda_\ell^m(\cos[\pi-\theta]) = (-1)^{\ell+m} \lambda_\ell^m(\cos\theta)\), which takes only a quick check to either negate a term or not.

All of these optimizations follow very closely the logic described in detail in Part I, so we’ll forego stepping through each of them individually. The two things to note are:

  1. Despite the sum in Eqn. \(\ref{eqn:fast_analysis}\) being written such that the implied outer iteration is over \(\ell\) and \(m\) and the sums over \(\theta\) and \(\phi\) are performed for each \((\ell,m)\), the nesting order of the coded loops does not change from the synthesis algorithm — we still iterate over pixels so that the FFT and per-ring \(\lambda_\ell^m\) optimizations apply, and instead we incrementally accumulate all \(a_{\ell m}\) simultaneously.

  2. Handling the case where \(m_\mathrm{max} \ge N_\phi\) requires aliasing coefficients much like with synthesis but in “reverse”. The alias_harmonic() function provided in Part I has been split into two functions — alias_index() and alias_coeffs() which together do the same operation, but splitting the implementations allows for calculating the appropriate aliased index before accessing the phases array.

The following function implements our fast spherical harmonic analysis for an ECP pixelization. Comparison to the matching synthesize_ecp function from the previous article should show that they are very similar, with only minor reordering and renaming of variables.

using FFTW

# Calculates aliased index `i′` given original index `i` and Fourier axis length 'len'.
# Also indicates whether conjugation or twice-real symmetries are needed.
@inline function alias_index(len::Int, i::Int)
    isconj, isnyq = false, false
    nyq = max(1, len ÷ 2)
    i < nyq && return (i, isconj, isnyq)
    i = mod(i, len)
    if i > nyq
        i = len - i
        isconj = true
    elseif i == 0 || (iseven(len) && i == nyq)
        isnyq = true
    end
    return (i, isconj, isnyq)
end
# Performs conjugation or twice-real transformation, if required, in mode aliasing.
@inline function alias_coeffs(coeffs::NTuple, isconj::Bool, isnyq::Bool)
    return ifelse(isnyq, complex.(2 .* real.(coeffs)),
                  ifelse(isconj, conj.(coeffs), coeffs))
end

function analyze_ecp(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}
    C = complex(R)
    ,  = size(map)
    nϕr =  ÷ 2 + 1 # real-symmetric FFT's Nyquist length (index)
    nθh = ( + 1) ÷ 2 # number of rings in northern hemisphere

    # pixel grid definition for ECP
    θr = crange(R(0.0), R(π), )
    ϕ₀ = R(π) / 
    ΔΩ = 2R(π)^2 / ( * )

    alms = fill(zero(C), lmax + 1, mmax + 1)
    Λ = zeros(R, lmax + 1, mmax + 1)

    @inbounds for j in 1:nθh
        j′ =  - j + 1 # southern ring index
        ,  = sincos(θr[j])
        λlm!(Λ, lmax, mmax, ) # λ_ℓ^m(cos θ) factors

        h₁ = rfft(map[j,:])  # phase factors for northern ring
        h₂ = rfft(map[j′,:]) # phase factors for southern ring

        sθΔΩ =  * ΔΩ
        for m in 0:mmax
            i, isconj, isnyq = alias_index(, m)
            a₁, a₂ = alias_coeffs((h₁[i+1], h₂[i+1]), isconj, isnyq)
            a₁, a₂ = (a₁, a₂) .* (sθΔΩ * cis(m * -ϕ₀))
            # Σ_{ℓ = m}^{lmax}
            for  in m:lmax
                c = isodd(+m) ? a₁ - a₂ : a₁ + a₂
                alms[+1,m+1] += c * Λ[+1,m+1]
            end
        end
    end
    return alms
end

As expected, simple timing of analyzing the map from Figure 2.1 shows that this implementation executes much quicker than the naive algorithm:

@time analyze_ecp(map, 100);
  0.034899 seconds (16.75 k allocations: 8.966 MiB)

and we observe that the optimized and reference functions give equivalent answers (within an acceptable difference in floating point values):

analyze_ecp(map, 100)  analyze_reference(map, 100)
true

Numerical accuracy of analyzing discretized maps

Let us return to the approximation made in Eqn. \(\ref{eqn:analysis_integral_sum}\) wherein we assumed one could replace the integral with a discrete summation. A trivial test is to vary the resolution of the map which is analyzed and see how the recovered \(a_{\ell m}\)s change. The example given above in Section II used a map which had been synthesized to a \(500 \times 1000\) pixel image, but we can easily reduce the resolution by a factor of \(5\). We do exactly this in Figure 4.1 below.

Similar to Figure 2.2 but shows a comparison between power spectra obtained from maps of two separate resolutions. The hi-res map is the same as above at \(500 \times 1000\) pixels, while the lo-res map is decreased by a factor of 5 to only \(100 \times 200\). We see here that the fractional error of analysis increases with decreasing resolution of the map.

It is clear that the lower resolution map results in much larger fractional differences than the original higher resolution map.

A simple test shows that having a complicated, varying continuous field discretized into pixels is not the problem, but rather that the approximation of the integral as a discrete sum itself is the underlying issue. Take the first spherical harmonic: \begin{align*} Y_{00}(\theta, \phi) &= \frac{1}{\sqrt{4\pi}} \end{align*} Since it is uniform everywhere, there is no place where a pixel can be “poorly located” to miss a change in the value on the sphere. Let us analyze this uniform map where we expect \(a_{00} = 1\) and everywhere else \(a_{\ell m} = 0\).

analyze_ecp(fill(1 / sqrt(4π), 500, 1000), 10, 3)
11×4 Array{Complex{Float64}, 2}:
    1.0000016449359603 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                   0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
 3.6782267450220785e-6 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                   0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
  4.934978353682631e-6 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                   0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
   5.93133121302037e-6 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                   0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
  6.783088214334931e-6 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                   0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
  7.539475913250632e-6 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im

What we see here is that not only is the \(a_{00}\) term greater than 1 (and by an amount where only the first 5 of 17 decimal digits are correct), we also see non-negligible response in the even \(\ell > 0\), \(m=0\) modes.

An interesting observation can be made if we now vary the resolution of the map — we can see that the unagreement between the ideal expectation and numerical results decrease as the number of pixels increases.

for npix in [100, 1_000, 10_000]
    alms_Y00 = analyze_ecp(fill(1/sqrt(4π), npix, 2npix), 10)
    println("npix = $(rpad(string(npix), 5)) : Δa00 = $(alms_Y00[1,1] - 1.0)")
end
npix = 100   : Δa00 = 4.112453549298678e-5 + 0.0im
npix = 1000  : Δa00 = 4.112336344785916e-7 + 0.0im
npix = 10000 : Δa00 = 4.112333806816082e-9 + 0.0im

By increasing the number of pixels by a factor of 100, the residual error in the value of \(a_{00}\) decreases in sympathy by a factor of 100.

While definitely not proof, it is at least suggestive that the discrepancy between numerical and analytical results can be explained by the discretization of the map and how it is being analyzed — extrapolating the trend suggests that infinite resolution would converge to the expected analytical result.

There turns out to be two common ways (at least that I’ve been introduced to) to improve the accuracy of the analysis step, and both are stated in the documentation for the analysis tool anafast from the HEALPix library:

Anafast permits two execution options which allow a significant improvement of accuracy of the approximate quadrature performed by this facility:

  • Improved analyses using either the provided ring weights, which correct the quadrature on iso-latitude rings, or pixel-based weights which improve the quadrature on every pixels, and/or
  • An iterative scheme using in succession several backward and forward harmonic transforms of the maps.

The second option provided is actually the simpler solution to implement, so we’ll start there.

Iterative quadrature

Resynthesizing the just-analyzed \(a_{\ell m}\)s and plotting the difference versus the original input is an illustrative exercise for why an iterative approach is motivated.

npix = 100
Y00_map = fill(1/sqrt(4π), npix, 2npix)
Y00_alm = analyze_ecp(Y00_map, 10)
Y00_resynth = synthesize_ecp(Y00_alm, npix, 2npix)
ΔY00 = Y00_map - Y00_resynth
Plot of the residual map \(\Delta Y_{00}\) which is derived from the difference between the analytic input \(Y_{00} = (4\pi)^{-1/2}\) and the resynthesis of an analyzed set \(\{a_{\ell m}\}\) calculated from a uniform discrete map of size \(100 \times 200\) pixels.

The residual difference map shows a pattern of zonal harmonics, which is qualitatively consistent with the observation above that analysis of the \(Y_{00}\) mode has non-zero even \(\ell\) (\(m=0\)) coefficients. The anafast documentation quoted above states that iterative reanalysis of residuals converges to the expected solution.

Let us do one iteration by hand as an illustrative example. Taking the residual map above (Figure 4.2):

  1. analyze the residual map back to a set of residual harmonic coefficients, \(\delta a_{\ell m}\).
  2. sum these with the previous set of coefficients: \(a_{\ell m} + \delta a_{\ell m} \rightarrow {}_{(2)}a_{\ell m}\)
  3. synthesize the new, iterated coefficients \({}_{(2)}a_{\ell m}\) back to a map
  4. compare the second iteration’s result with the original input map

The following code does this iteration once, and Figure 4.3 shows the result.

ΔY00_alm = analyze_ecp(ΔY00, 10)
Y00_resynth2 = synthesize_ecp(Y00_alm + ΔY00_alm, npix, 2npix)
ΔY00_iter2 = Y00_map - Y00_resynth2
The residual remaining after one round of iteration — as described in the text above — beyond Figure 4.2 of the analysis of the uniform map \(Y_{00} = (4\pi)^{-1/2}\). The scale is the same as in Figure 4.2, so the uniform gray shows that the residual is significantly closer to zero than the analysis without any iteration.

We see that a single iteration has decreased the remaining residual closer to zero than can be easily seen with the fixed color scale used here. (The color scales in both Figures 4.2 and 4.3 are the same to be a fair comparison.)

Rather than doing these iterations by hand, we can easily write a wrapper function which does the iteration automatically, continuing until either (a) the maximum absolute residual is smaller than some tolerance, or (b) a maximum number of iterations have been completed.

function analyze_ecp_iter(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax;
                          atol = sqrt(eps(maximum(abs, map))),
                          maxiter::Integer = 100) where {R<:Real}
    n, m = size(map)
    alms = analyze_ecp(map, lmax, mmax)
    for _ in 2:Int(maxiter)  # bail if maxiter reached
        map′ = synthesize_ecp(alms, n, m)
        δmap = map - map′
        # check for convergence
        if maximum(abs, δmap) < atol
            break
        end
        alms .+= analyze_ecp(δmap, lmax, mmax)
    end
    return alms
end

With this function, we can again test for how closely analysis of the \(Y_{00}\) mode results in \(a_{00} = 1\),

for npix in [100, 1_000, 10_000]
    alms_Y00 = analyze_ecp_iter(fill(1/sqrt(4π), npix, 2npix), 10)
    println("npix = $(rpad(string(npix), 5)) : Δa00 = $(alms_Y00[1,1] - 1.0)")
end
npix = 100   : Δa00 = 3.063957976223719e-10 + 0.0im
npix = 1000  : Δa00 = -1.1162182289581324e-11 + 0.0im
npix = 10000 : Δa00 = -9.992007221626409e-16 + 0.0im

and we see that the deviation from expectation has decreased by many orders of magnitude. In the last case, we’ve actually achieved an error of just 4.5 ulps, and doing so required nothing more sophisticated than an additional loop.

Weighted quadrature

The other method which is mentioned by anafast is to make use of “ring weights” or “pixel weights” to reduce the integration error.

The generic premise is that the integral can be better approximated by a summation if we include some weight factors \(w_i\) at each integration point \(x_i\) so that \begin{align} \int f(x)\,dx \approx \sum_i w_i f(x_i) \end{align} in such a way that the summation is a better approximation to the integral than for the fixed “weight” \(w_i = \Delta x\).

The aforementioned HEALPix tools provide the quadrature weights for their pixelization scheme, but deriving such generalized weights is beyond the scope of this article. Instead, we can take the approach of using the well-known Gauss-Legendre quadrature rule which requires that the locations of the pixel centers be changed to agree with the given weights. (This may not be a reasonable trade-off in a real application, but it makes the discussion and implementation easier here.)

For the sake of brevity, I’ll skip over describing the necessary modifications to the synthesis/analysis functions, but hopefully the changes hidden within the “Show GL-node code” button to the below-right are straight-forward to follow. The important details to note are that the new functions synthesize_gl and analyze_gl no longer produce/consume maps on an ECP grid, but rather they operate on a grid where

  1. the azimuthual pixel centers are unchanged — uniformly spaced grid from \(0\) to \(2\pi\)
  2. the colatitudinal pixel centers are moved to the nodes of the Gauss-Legendre quadrature rule

and the corresponding quadrature weight is applied for each ring in the map.

The calculation of the nodes and weights for the Gauss-Legendre quadrature used in the synthesis and analysis functions below is performed by FastGaussQuadrature.jl.

using FastGaussQuadrature: gausslegendre

function synthesize_gl(alms::Matrix{C}, , ) where {C<:Complex}
    lmax, mmax = size(alms) .- 1
    z, _ = gausslegendre() # cos(θ) and corresponding weights
    ϕ₀ = π /  # == first(crange(0.0, 2.0π, nϕ))
    nϕr =  ÷ 2 + 1 # length real-only half of FFT axis; see `rfft()`
    nθh = ( + 1) ÷ 2 # number of rings in northern hemisphere

    R = real(C)
    Λ = Matrix{R}(undef, size(alms)...)
    g₁ = zeros(C, nϕr) # phase factors for northern ring
    g₂ = zeros(C, nϕr) # phase factors for southern ring
    map = zeros(R, , )

    for y in 1:nθh
        y′ =  - y + 1 # southern ring index
        λlm!(Λ, lmax, mmax, -z[y]) # λ_ℓ^m(cos θ) factors; N.B. negate z
                                   # because θ and z increase in opposite directions
        for m in 0:mmax
            accn, accs = zero(C), zero(C)
            # Σ_{ℓ = m}^{lmax}
            for  in m:lmax
                term = alms[+1,m+1] * Λ[+1,m+1]
                accn += term
                accs += isodd( + m) ? -term : term
            end
            accn, accs = (accn, accs) .* cis(m * ϕ₀) # ring offset rotation
            i, isconj, isnyq = alias_index(, m)
            accn, accs = alias_coeffs((accn, accs), isconj, isnyq)
            g₁[i+1] += accn
            g₂[i+1] += accs
        end
        map[y, :] = brfft(g₁, )
        fill!(g₁, zero(C))
        map[y′, :] = brfft(g₂, )
        fill!(g₂, zero(C))
    end
    return map
end

function analyze_gl(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}
    C = complex(R)
    ,  = size(map)
    nϕr =  ÷ 2 + 1 # real-symmetric FFT's Nyquist length (index)
    nθh = ( + 1) ÷ 2 # number of rings in northern hemisphere

    # pixel grid definition for gauss-legendre grid
    z, w = gausslegendre() # cos(θ) and corresponding weights
    ϕ₀ = R(π) /  # == first(crange(0.0, 2.0π, nϕ))
    Δϕ = 2R(π) / 

    alms = fill(zero(C), lmax + 1, mmax + 1)
    Λ = zeros(R, lmax + 1, mmax + 1)
    f = zeros(R, )

    @inbounds for j in 1:nθh
        j′ =  - j + 1 # southern ring index
        λlm!(Λ, lmax, mmax, -z[j]) # λ_ℓ^m(cos θ) factors; N.B. negate z
                                   # because θ and z increase in opposite directions

        h₁ = rfft(map[j,:])  # phase factors for northern ring
        h₂ = rfft(map[j′,:]) # phase factors for southern ring

        wΔϕ = w[j] * Δϕ  # N.B. no sin(θ) factor due to change of variables
                         #        z = cos(θ) -> dz = sin(θ) dθ
                         #      and w implicitly accounts for that factor.
        for m in 0:mmax
            i, isconj, isnyq = alias_index(, m)
            a₁, a₂ = alias_coeffs((h₁[i+1], h₂[i+1]), isconj, isnyq)
            a₁, a₂ = (a₁, a₂) .* (wΔϕ * cis(m * -ϕ₀))
            # Σ_{ℓ = m}^{lmax}
            for  in m:lmax
                c = isodd(+m) ? a₁ - a₂ : a₁ + a₂
                alms[+1,m+1] += c * Λ[+1,m+1]
            end
        end
    end
    return alms
end

Repeating the test of analyzing the \(Y_{00}\) uniform mode across the sphere:

analyze_gl(fill(1 / sqrt(4π), 500, 1000), 10, 3)
11×4 Array{Complex{Float64}, 2}:
      0.9999999999999997 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                     0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
 -1.4224732503009818e-16 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                     0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
   -7.37257477290143e-17 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                     0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
  1.1709383462843448e-16 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                     0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
 -1.0408340855860843e-17 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
                     0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im
  1.2663481374630692e-16 + 0.0im  0.0 + 0.0im  0.0 + 0.0im  0.0 + 0.0im

we see that unlike the (uniterated) ECP grid, we achieve very high numerical accuracy in recovering the expected spherical harmonic coefficients. In fact, lowering the resolution of the grid has negligible impact on how accurately the value is recovered:

for npix in [100, 1_000, 10_000]
    alms_Y00 = analyze_gl(fill(1/sqrt(4π), npix, 2npix), 10)
    println("npix = $(rpad(string(npix), 5)) : Δa00 = $(alms_Y00[1,1] - 1.0)")
end
npix = 100   : Δa00 = 2.220446049250313e-16 + 0.0im
npix = 1000  : Δa00 = -4.440892098500626e-16 + 0.0im
npix = 10000 : Δa00 = 2.220446049250313e-16 + 0.0im

which clearly emphasizes how important the issue of quadrature weights can be when high accuracy results are needed.

Conclusions

This article completes the discussion — started in Part I — in demonstrating the reverse operation of analyzing discrete maps on the sphere to determine the spherical harmonic coefficients. On top of the symmetry and aliasing issues previously explored, we ecountered a new concern related to the obligatory step of approximating an integral with sums.

We saw that a simple solution is to just iterate a pair of synthesis and reanalysis of the residual until a desired level of convergence is achieved, which is an attractive option because it requires no further specialized knowledge. The obvious disadvantage is that a large number of back-and-forward spherical harmonic transforms requires a significant amount of extra computational resources.

The alternative which we briefly introduced is to instead make use of quadrature rules to carefully construct a situation where modifying the pixelization and combining with a set of numerical weight factors greatly increases the numerical precision of the analysis without requiring any iteration. The cost here, though, is that (at least in the example case described here) the pixel grid is no longer our choice but rather has to follow the rules of the quadrature scheme.