Pre-normalizing Legendre Polynomials

In the previous part we computed Legendre polynomials via a few recurrence relations, but the polynomial values grow rapidly with degree and quickly overflow the range of standard finite floating point numbers. In this posting we explore baking in a normalization factor into the recurrence relation — such as used when the Legendre polynomials are used to calculate the spherical harmonics — that eliminates the overflow.

Introduction¶

At the end of the previous article we saw that the growth of the Legendre polynomials polynomials rapidly leads to numerical overflow. For common uses such as calculating the spherical harmonics in CMB analysis, we require calculating the polynomials to much higher degrees than where the overflows occur, so our goal in this article is to explain how to mitigate this problem.

As a reminder, recall that we calculate the associated Legendre polynomials by evaluating the following three recurrence relations — starting from an appropriate initial condition — for any integer degree $\ell$ and order $m$.

\begin{align} P_\ell^\ell(x) &= -(2\ell-1)\sqrt{1-x^2} P_{\ell-1}^{\ell-1}(x) \label{eqn:std_rr_1term_lm} \\ P_\ell^{\ell-1}(x) &= (2\ell-1) x P_{\ell-1}^{\ell-1}(x) \label{eqn:std_rr_1term_l} \\ P_\ell^m(x) &= \frac{2\ell-1}{\ell - m} x P_{\ell-1}^m(x) - \frac{\ell+m-1}{\ell - m} P_{\ell-2}^m(x) \label{eqn:std_rr_2term} \end{align}

Note that compared to Eqs. 1–3 in the previous article, we have shifted the indexing so that $\ell+1 \rightarrow \ell$ and the left-hand side is for degree $\ell$ (not $\ell+1$). (This will prove to be more convenient shortly when we start generalizing these equations.)

This article will be structured in two primary parts. In the first part we’ll generalize the recurrence relations to prepare for supporting a wider range of normalizations. Then in the second part we’ll show how to pre-normalize the Legendre functions for use in the spherical harmonics, avoiding the overflow which would otherwise prohibit complete analysis of CMB skies.

Generalized recurrence relations¶

We’ll rewrite the recurrence relations for a generalized form of the Legendre polynomials $p_\ell^m(x)$ which include some arbitrary normalization factor in terms of four undetermined coefficients $\mu_\ell$, $\nu_\ell$, $\alpha_\ell^m$, and $\beta_\ell^m$:

\begin{align} p_\ell^\ell(x) &= -\mu_\ell \sqrt{1-x^2} p_{\ell-1}^{\ell-1}(x) \label{eqn:cus_rr_1term_lm} \\ p_\ell^{\ell-1}(x) &= \nu_\ell x p_{\ell-1}^{\ell-1}(x) \label{eqn:cus_rr_1term_l} \\ p_\ell^m(x) &= \alpha_\ell^m x p_{\ell-1}^m(x) - \beta_\ell^m p_{\ell-2}^m(x) \label{eqn:cus_rr_2term} \end{align}

For the standard associated Legendre polynomials where $P_0^0(x) = 1$ (which will from here on out call the unit-normalized Legendre polynomials), we can simply read off the coefficients from the original form of the recurrence relations and identify each of the coefficients as:

\begin{align}\begin{aligned} \mu_\ell &= 2\ell - 1 & \alpha_\ell^m &= \frac{2\ell - 1}{\ell - m} \\ \nu_\ell &= 2\ell - 1 & \beta_\ell^m &= \frac{\ell + m - 1}{\ell - m} \end{aligned}\end{align}

To prepare for a fully generalized implementation, we’ll encode the calculation of each coefficient term and the initial condition into functions which are dispatched on a trait type¹. At this point the addition of a trait type might seem unnecessary, but the use of traits will let us use a single implementation to calculate any arbitrary normalization².

# Unit-normalization trait type
struct UnitNorm end

# Initial condition and recurrence relation coefficients for unit normalization
initcond(::UnitNorm) = 1.0
coeff_μ(::UnitNorm, l) = 2l - 1
coeff_ν(::UnitNorm, l) = 2l - 1
coeff_α(::UnitNorm, l, m) = (2l - 1) / (l - m)
coeff_β(::UnitNorm, l, m) = (l + m - 1) / (l - m)

Now we update the implementation of legendre() that we wrote in the previous article to instead call these helper functions instead of having the initial condition and recurrence coefficients written inline. Note that at the moment we hard-code calls to use the unit normalization (such as in initcond(UnitNorm())), so the implementation is not yet generic, but this gets us a step closer.

function legendre(ℓ, m, x)
    P = convert(typeof(x), initcond(UnitNorm()))
    ℓ == 0 && m == 0 && return P

    # Use recurrence relation 1 to increase to the target m
    y = sqrt(one(x) - x^2)
    for m′ in 1:m
        μ = coeff_μ(UnitNorm(), m′)
        P = -μ * y * P
    end
    ℓ == m && return P

    # Setup for 2-term recurrence, using recurrence relation 2 to set
    # Pmin1
    ν = coeff_ν(UnitNorm(), m + 1)
    Pmin2 = P
    Pmin1 = ν * x * Pmin2
    ℓ == m+1 && return Pmin1

    # Use recurrence relation 3
    for ℓ′ in m+2:ℓ
        α = coeff_α(UnitNorm(), ℓ′, m)
        β = coeff_β(UnitNorm(), ℓ′, m)
        P = α * x * Pmin1 - β * Pmin2
        Pmin2, Pmin1 = Pmin1, P
    end

    return P
end

Comparison of legendre() from the previous posting versus the code cell above.³

 function legendre(ℓ, m, x)
-    P = one(x)
+    P = convert(typeof(x), initcond(UnitNorm())
     ℓ == 0 && m == 0 && return P
 
     # Use recurrence relation 1 to increase to the target m
     y = sqrt(one(x) - x^2)
     for m′ in 1:m
-        P = -(2m′ - 1) * y * P
+        μ = coeff_μ(UnitNorm(), m′)
+        P = -μ * y * P
     end
     ℓ == m && return P
 
     # Setup for 2-term recurrence, using recurrence relation 2 to set
     # Pmin1
+    ν = coeff_ν(UnitNorm(), m + 1)
     Pmin2 = P
-    Pmin1 = (2*(m+1) - 1) * x * Pmin2
+    Pmin1 = ν * x * Pmin2
     ℓ == m+1 && return Pmin1
 
     # Use recurrence relation 3
     for ℓ′ in m+2:ℓ
-        P = ((2ℓ′ - 1) * x * Pmin1 - (ℓ′ + m - 1) * Pmin2) / (ℓ′ - m)
+        α = coeff_α(UnitNorm(), ℓ′, m)
+        β = coeff_β(UnitNorm(), ℓ′, m)
+        P = α * x * Pmin1 - β * Pmin2
         Pmin2, Pmin1 = Pmin1, P
     end
 
     return P
 end

No functionality has yet actually changed for an end-user of the function. We can exactly reproduce the figure from last time showing the $m = 1$ associated Legendre polynomials:

using PyPlot

x = range(-1.0, 1.0, length=1000)
for m in 1:5
    plot(x, legendre.(m, 1, x), label="\$P_{$m}^1\$")
end
axhline.([-1, 0, 1], color="gray", linewidth=0.5)
axvline.([-1, 0, 1], color="gray", linewidth=0.5)
legend(loc = "lower right")
title("Associated Legendre Polynomials (\$m = 1\$)")

We’re now ready to tackle the issue of the Legendre polynomials growing in magnitude with increasing degree and order by considering how the spherical harmonics are normalized.

Legendre Polynomials in the Spherical Harmonics¶

The spherical harmonics are critical to CMB analysis, so accurately and efficiently calculating them is worth the effort. In part 1 we defined the spherical harmonics as

\begin{align} Y_{\ell m}(\theta,\phi) = \sqrt{\frac{2\ell+1}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}} P_\ell^m(\cos\theta) e^{im\phi} \end{align}

which composes the Legendre polynomials (as a function of colatitude $\theta$) with a complex exponential (acting on the azimuth $\phi$) and a normalization constant. We then saw in the previous part that the Legendre polynomials result in values greater than what is representable in standard floating point types at very modest combinations of degree and order (around $\ell,m \sim 150$ for 64-bit floating point numbers).

legendre(152, 150, 0.2)

Inf

We discussed last time that the normalization counteracts growth rate of the polynomials, but it’s not sufficient to calculate the polynomials and normalization independently. In fact, the normalization underflows to identically zero when we know that mathematically it’s non-zero (just very small):

Nlm(l, m; type=Float64) = sqrt((2l + 1) / 4π / reduce(*, l-m+1:l+m, init=one(type)))

Nlm(152, 150)

0.0

A big feature of Julia is that it specializes your code depending on the types of the input arguments, and that allows one workaround for under/overflowing to be to upgrade to using extended-precision BigFloat values instead of standard Float64 types. With BigFloat, we can calculate both the polynomial value and normalization and then convert back to Float64:

using Printf
n = Nlm(152, 150, type=BigFloat)
p = legendre(152, 150, big(0.2))
λ = Float64(n * p)
@printf """
    N_152^150 ≈ %0.9g
    P_152^150 ≈ %0.9g
    λ_152^150 ≈ %0.9g
    """ n p λ

N_152^150 ≈ 1.32090549e-309
P_152^150 ≈ 2.94031628e+308
λ_152^150 ≈ 0.388387991

Using BigFloat comes at a heavy cost, though. Compare the computation time for a standard 64-bit machine float versus calculating on BigFloat values:

using BenchmarkTools

print(" Float64:"); @btime legendre(152, 150, 0.2);
print("BigFloat:"); @btime legendre(152, 150, big(0.2));

 Float64:  146.484 ns (0 allocations: 0 bytes)
BigFloat:  20.763 μs (624 allocations: 34.13 KiB)

The calculation using BigFloat values takes over two orders of magnitude longer to calculate than the equivalent Float64 value. Doing similar calculations for thousands of degrees and orders, for potentially thousands more input values, means the extra cost is unacceptable if there’s any other option.

Spherical Harmonic Normalization¶

Since it’s been promised already, the obvious answer is that we can instead bake-in the normalization into the recurrence relations so that we calculate the final value directly.

The first step is to note that we can define the spherical harmonic pre-normalized Legendre polynomials $\lambda_\ell^m(x)$ as

\begin{align} \lambda_\ell^m(x) \equiv N_\ell^m P_\ell^m(x) \quad\text{where}\quad N_\ell^m \equiv \sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}} \end{align}

by combining the normalization and Legendre polynomials into a single function so that the spherical harmonics are rewritten simply as

\begin{align} Y_{\ell m}(\theta,\phi) = \lambda_\ell^m(\cos\theta) e^{im\phi} \end{align}

Our remaining task is to determine the expressions for the coefficients in Eqs. $\ref{eqn:cus_rr_1term_lm}$–$\ref{eqn:cus_rr_2term}$ that produce the pre-normalized functions $\lambda_\ell^m(x)$. We do that by multiplying the original [unit-normlized] recurrences (Eqs. $\ref{eqn:std_rr_1term_lm}$–$\ref{eqn:std_rr_2term}$) through by $N_\ell^m$. The left-hand side will by definition be $\lambda_\ell^m$, and on the right-hand side we rearrange terms until we can similarly identify the $\lambda_{\ell’}^{m’}$ terms.

For example, $\alpha_\ell^m$ and $\beta_\ell^m$ are determined by starting with Eq. $\ref{eqn:cus_rr_2term}$. Multiplying through by $N_\ell^m$ and identifying $\lambda_\ell^m$ on the left-hand side, that leaves us with

\begin{align} \begin{split} \lambda_\ell^m &= \frac{2\ell-1}{\ell-m} x \sqrt{\frac{2\ell+1}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}} P_{\ell-1}^m(x) - \frac{\ell+m-1}{\ell-m} \sqrt{\frac{2\ell+1}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}} P_{\ell-2}^m(x) \end{split} \end{align}

Through judicious use of algebra, the terms on the right-hand side can be manipulated to gather terms of the form $N_{\ell-1}^m P_{\ell-1}^m(x)$ and $N_{\ell-2}^m P_{\ell-2}^m(x)$, leaving us with

\begin{align} \lambda_\ell^m &= \sqrt{\frac{2\ell+1}{2\ell-3} \frac{4(\ell-1)^2 - 1}{\ell^2 - m^2}} x \lambda_{\ell-1}^m(x) - \sqrt{\frac{2\ell+1}{2\ell-3} \frac{(\ell-1)^2 - m^2}{\ell^2 - m^2}} \lambda_{\ell-2}^m(x) \end{align}

We identify each of the two square root terms as $\alpha_\ell^m$ and $\beta_\ell^m$ since they are the cofficients appropriate for generating $\lambda_\ell^m(x)$.

Proceeding similarly with the other two recurrence relations, we obtain:

\begin{align}\begin{aligned} \mu_\ell &= \sqrt{1 + \frac{1}{2\ell}} & \alpha_\ell^m &= \sqrt{\frac{2\ell+1}{2\ell-3} \frac{4(\ell-1)^2 - 1}{\ell^2 - m^2}} \\ \nu_\ell &= \sqrt{1 + 2\ell} & \beta_\ell^m &= \sqrt{\frac{2\ell+1}{2\ell-3} \frac{(\ell-1)^2 - m^2}{\ell^2 - m^2}} \end{aligned}\end{align}

Finally, we also require the intial condition $\lambda_0^0(x)$. This too is straight forward given the definition:

\begin{align} \lambda_0^0(x) &= N_0^0 P_0^0(x) = \sqrt{\frac{1}{4\pi}} \end{align}

Computing this normalization is accomplished by definining a new trait type SphereNorm and overloading the initial condition and coefficient functions defined above, dispatching on SphereNorm.

struct SphereNorm end

# Initial condition and recurrence relation coefficients for unit normalization
initcond(::SphereNorm) = 1.0 / sqrt(4π)
coeff_μ(::SphereNorm, l) = sqrt(1 + 1/2l)
coeff_ν(::SphereNorm, l) = sqrt(1 + 2l)
coeff_α(::SphereNorm, l, m) = sqrt((2l+1)/(2l-3) * (4*(l-1)^2 - 1)/(l^2 - m^2))
coeff_β(::SphereNorm, l, m) = sqrt((2l+1)/(2l-3) * ((l-1)^2 - m^2)/(l^2 - m^2))

To make use of the new normalization, we now simply update legendre() to accept the normalization as an argument and replace each use of UnitNorm() with the normalization argument.

function legendre(norm, ℓ, m, x)
    P = convert(typeof(x), initcond(norm))
    ℓ == 0 && m == 0 && return P

    # Use recurrence relation 1 to increase to the target m
    y = sqrt(one(x) - x^2)
    for m′ in 1:m
        μ = coeff_μ(norm, m′)
        P = -μ * y * P
    end
    ℓ == m && return P

    # Setup for 2-term recurrence, using recurrence relation 2 to set
    # Pmin1
    ν = coeff_ν(norm, m + 1)
    Pmin2 = P
    Pmin1 = ν * x * Pmin2
    ℓ == m+1 && return Pmin1

    # Use recurrence relation 3
    for ℓ′ in m+2:ℓ
        α = coeff_α(norm, ℓ′, m)
        β = coeff_β(norm, ℓ′, m)
        P = α * x * Pmin1 - β * Pmin2
        Pmin2, Pmin1 = Pmin1, P
    end

    return P
end

With the normalization now baked into the recurrence, we have no problem calculating the appropriately-normalized polynomial value for $\lambda_{152}^{150}(0.2)$.

legendre(SphereNorm(), 152, 150, 0.2)

0.3883879907461426

It could be seen as unfortunate that the call syntax requires we explicitly state which normalization is to used, but it’s trivial to define wrapper functions that are named for each normalization style. In accordance with the mathematical notation of this article, we could choose to define Plm and λlm for the unit and spherical harmonic normalizations, respectively.

Plm(l, m, x) = legendre(UnitNorm(), l, m, x)
λlm(l, m, x) = legendre(SphereNorm(), l, m, x)

Using the pre-normalization, we can now update the plot from the previous posting which had missing values due to overflow in the middle and show the two functions entirely:

plot(x, λlm.(151, 150, x), label="\$λ_{151}^{150}\$")
plot(x, λlm.(152, 150, x), label="\$λ_{152}^{150}\$")
legend()
title("Spherical Harmonic normalized Legendre polynomials")

Conclusions¶

The trait-based dispatch mechanism is one of the core features of the Legendre.jl package, and we could undoubtedly use this simple implementation in real use cases. To achieve highly efficient and high-accuracy results, though, the package implements several other optimization and numerical tricks which we will introduce and discuss in future postings. We’ll explore some of these enhancements as well, using several goal calculations to demonstrate existing deficiencies and how I have resolved them.

Traits in Julia are sometimes also called “Holy Traits” in recognition of a core Julia contributor, Tim Holy, who early-on recognized the power of Julia’s multiple dispatch combined with singleton instance types as a way to make extensible code and type hierarchies. Trait types are now used extensively throughout base Julia, such as in defining the iterator interface (e.g. HasLength(), HasEltype(), etc). ↩︎
The implementation in this article is simplified from the version actually used in Legendre.jl. For instance, we’ve skipped using an abstract supertype to unify all normalizations under one type hierarchy. Furthermore, the interface functions have all been simplified with the assumption that we’re only interested in Float64 results; supporting type-stable results over all numeric types requires passing additional contextual information to the calculations. ↩︎
The diff is computed against a version of legendre() which has had the indexing updated for the $\ell+1 \rightarrow \ell$ change, so the before state does not exactly correspond to the legendre() function shown in the previous posting. This has been done to highlight the trait-based interface related changes. ↩︎