The Associated Legendre Polynomials are an important set of functions in cosmic microwave background (CMB) research. In this first part of an upcoming series, I motivate the need for a new high-performance Julia package — Legendre.jl — which I am writing to recreate a key algorithm in my thesis research.

Introduction¶

My graduate work with the Bicep/Keck Array Collaboration researched the cosmic microwave background (CMB), the afterglow of the Big Bang. We were aimed at finding evidence for a hypothesized, very brief period of exponential expansion of the Universe called Inflation which would nicely explain modern cosmological observations. The evidence would come in the form of observing so-called “B-mode” patterns of polarization in the CMB. The challenge, though, is that B-modes are much dimmer than a different kind of polarization pattern — E-modes — that dominate any attempt at extracting the B-mode signal if no mitigation is made during analysis.

A core component of my thesis¹ revolved around work I did in that mitigation effort, making use of a matrix-based algorithm² for “purifying” the maps we constructed of the CMB before further processing them. One input to calculating the purification operator is constructing a “pixel-pixel covariance matrix” which encodes how points on the sky are correlated with one another, and calculating this covariance matrix is a numerically expensive operation that involves evaluating the Associated Legendre Polynomials many times.

In this article, I will introduce the Associated Legendre Polynomials and motivate why I have decided to write a new implementation in Julia, available in the Legendre.jl package. Later articles in this series will discuss in more detail specific aspects of the implementation I have discovered along the way as well as explore various use cases for a generic Legendre calculation package.

Spherical Harmonics¶

Before jumping to the pixel-pixel covariance calculation, let us beging with a description every physicist (and probably most mathematicians) will already be familiar with — the spherical harmonics. For a review of deriving the spherical harmonics in detail, to physicists I’d suggest standard texts such as Griffiths’ Introduction to Quantum Mechanics³ or Jackson’s Classical Electrodynamics⁴ which they probably already have on hand, but any text that deals with the solution of second-order differential equations in spherical coordinates should do.

The spherical harmonics arise from solving Laplace’s equation \begin{align} \nabla^2 \psi = 0 \label{eqn:wave_eqn} \end{align} in spherical coordinates. The equation is separable into a radial component $R(r)$ and an angular part $Y(\theta,\phi)$ such that the total solution is $\psi(r,\theta,\phi) \equiv R(r)Y(\theta,\phi)$. We’ll ignore the radial component since the CMB exists on the surface of the sphere and continue with only the angular part.

When arranged into the standard form (following the physicists' convention⁵), the further separated differential equations are

where $x = \cos\theta$ and the indices $\ell$ and $m$ are separation constants subject to the constraints that they the are both integers, $\ell \ge 0$, and $|m| \le \ell$. The solution to Eq. $\ref{eqn:diffeq_phi}$ are the complex exponentials, $\Phi_m(\phi) = e^{im\phi}$, and the Associated Legendre Polynomials (of the First Kind) $P_\ell^m(x)$ of degree $\ell$ and order $m$ are a solution to Eq. $\ref{eqn:diffeq_theta}$. Each are complete orthogonal sets in the azimuthal and polar coordinates, respectively, so together they form a complete orthogonal set on the surface of the sphere. The normalized spherical harmonics are therefore \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}

As a complete orthogonal set of functions, the spherical harmonics can be used to decompose functions on the sphere, much like how the complex exponentials are the basis functions of the Fourier Transform. For some function $f(\theta,\phi)$ on the sky (such as the CMB temperature anisotropies), the function can be decomposed into the set of harmonic coefficients $a_{\ell m}$ according to \begin{align} f(\theta,\phi) &= \sum_{\ell=0}^\infty \sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m}(\theta,\phi) & \text{where } a_{\ell m} &\equiv \int_0^{2\pi} \int_0^\pi f(\theta,\phi) Y_{\ell m}^*(\theta,\phi) ,\sin\theta ,d\theta,d\phi \end{align}

Both the forward and reverse transformation require calculating the Associated Legendre Polynomials for many orders, degrees, and arguments, so this common use motivates the need for an efficient means of calculating the functions numerically.

My Motivation¶

While very much critical for CMB analysis, the spherical harmonic transform was not the [direct] reason for my interest in calculating the Associated Legendre Polynomials. As I mentioned in the introduction, I was actually concerned with calculating a covariance matrix for use in a purification operation. A later article will discuss the details of the calculation, but one of the steps involes calculating the following four weight functions, which are themselves functions of the Legendre polynomials:

As a (very underestimated) lower bound on the computation time required, let us assume that each of the four equations can be completed in a single cycle of a processor. The argument $z_{ij}$ is the angular distance between pairs of points, and for my case, there were $\mathcal{O}\left(10^5\right)$ pixel locations. That means we must evaluate the expression for about $10^{10}$ values of $z_{ij}$. The expressions are also evaluated for a span of $\ell$ from $0$ to some $\ell_{\mathrm{max}}$, so assuming $\ell_{\mathrm{max}} = \mathcal{O}(10^3)$, that brings the total to (at minimum) $10^{10} \times 10^3 \times 4 = 4 \times 10^{13}$ cycles.

Ignoring the fact that each of the above equations can clearly not be computed in a single cycle, that still translates to about 5.5 hours for a single $2\text{ GHz}$ processor. Add in the extra costs of loading values from RAM into the CPU cache, requiring time to do all of the sub-expression additions/multiplications shown above, and evalutating the Associated Legendre Polynomials themselves (which we will later see also require a large number of sub-calculations), and the time could easily balloon by a couple of orders of magnitude if not implemented efficiently.

When I joined the collaboration, a working implementation of the entire purification process had recently been added to our processing pipeline, but it — like the rest of the analysis pipeline — was implemented in Matlab. In trying to take the what was initially designed as proof-of-concept code to a more optimized and generalized solution, I quickly concluded that the inability to control Matlab’s memory management was a bottleneck in further performance improvements. I ultimately replaced the core calculations with a C implementation that was accessible via MEX functions, but that translation is what initially put me on the path of learning how to calculate the Legendre polynomials.

The mathematical and algorithmic foundations of the Bicep/Keck Array purification process have all been previously published, but the only existing implementation exists as part of a private Matlab pipeline. While working day-to-day in Matlab for my graduate research, I have personally gravitated toward Julia as my programming language of choice for scientific computing. To that end, my goal is to provide the public with a Julia implementation of the Bicep/Keck Array purification process which is completely independent of our private Matlab implementation.

The beginning of that is writing a Julia implementation of calculating the Associated Legendre Polynomials, which I have done in Legendre.jl. The following installments in this article series will be focused on sharing the details being the package, especially pointing out areas where there were (at least to me) non-obvious concerns that I have had to work around. I’ll also take the opportunity at times to dicuss auxiliary uses of Legendre.jl that showcase it’s uses outside of calculating the pixel-pixel covariance matrices.