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^{1} revolved around work I did in that
mitigation effort, making use of a matrix-based algorithm^{2} 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*^{3} or
Jackson’s *Classical Electrodynamics*^{4} 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^{5}), the further separated differential equations are
`\begin{align} \frac{d^2}{d\phi^2}\Phi_m(\phi) &+ m^2\Phi_m(\phi) = 0 \label{eqn:diffeq_phi} \\ \left(1-x^2\right) \frac{d^2}{dx^2} P_\ell^m(x) &- 2x \frac{d}{dx}P_\ell^m(x) + \left[ \ell(\ell+1) - \frac{m^2}{1-x^2} \right] P_\ell^m(x) = 0 \label{eqn:diffeq_theta} \end{align}`

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¶

`\(\newcommand{\expv}[1]{\left\langle #1\right\rangle} \newcommand{\covF}[1]{F_\ell^{#1}(z_{ij})} \newcommand{\covTT}{\expv{T_i T_j}} \newcommand{\covTQ}{\expv{T_i Q_j}} \newcommand{\covTU}{\expv{T_i U_j}} \newcommand{\covQQ}{\expv{Q_i Q_j}} \newcommand{\covQU}{\expv{Q_i U_j}} \newcommand{\covUU}{\expv{U_i U_j}} \newcommand{\cij}{c_{ij}} \newcommand{\sij}{s_{ij}} \newcommand{\cji}{c_{ji}} \newcommand{\sji}{s_{ji}} \newcommand{\bigO}{\mathcal{O}}\)`

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:
`\begin{align*} \begin{aligned} \covF{00} &= \eta_\ell P_\ell(z_{ij}) \label{eqn:theorycov:F00} \\ \covF{10} &= \alpha_\ell \left[ \frac{z_{ij}}{1-z_{ij}^2} P_{\ell-1}(z_{ij}) - \left( \frac{1}{1-z_{ij}^2} + \frac{\ell-1}{2} \right)P_\ell(z_{ij}) \right] \label{eqn:theorycov:F10} \\ \covF{12} &= \beta_\ell \left[ \frac{\ell+2}{1-z_{ij}^2} z_{ij} P^2_{\ell-1}(z_{ij}) - \left( \frac{\ell-4}{1-z_{ij}^2} + \frac{\ell(\ell-1)}{2} \right) P^2_\ell(z_{ij}) \right] \label{eqn:theorycov:F12} \\ \covF{22} &= 2\beta_\ell \left[ \frac{\ell+2}{1-z_{ij}^2} P^2_{\ell-1}(z_{ij}) - \frac{\ell-1}{1-z_{ij}^2} z_{ij} P^2_\ell(z_{ij}) \right] \label{eqn:theorycov:F22} \end{aligned} \quad\text{where}\quad \begin{aligned} \eta_\ell &\equiv \frac{2\ell+1}{4\pi} \\ \alpha_\ell &\equiv \frac{2\eta_\ell \ell} {\sqrt{(\ell-1)\ell(\ell+1)(\ell+2)}} \\ \beta_\ell &\equiv \frac{2\eta_\ell} {(\ell-1)\ell(\ell+1)(\ell+2)} \end{aligned} \end{align*}`

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 `\(\bigO\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}} = \bigO(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\,\mathrm{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.

J. Willmert. “Constraining Inflationary

*B*-modes with the Bicep/*Keck Array*Telescopes”. PhD thesis. University of Minnesota, Nov 2019. url: http://hdl.handle.net/11299/211821 ↩︎The Bicep/

*Keck Array*Collaborations. “Bicep/*Keck Array*VII: Matrix Based*E*/*B*Separation Applied to Bicep2 and the*Keck Array*”. In: Astrophysical Journal 825, 66 (July 2016), p. 66. doi: 10.3847/0004-637X/825/1/66 arXiv: 1603.05976 [astro-ph.CO] ↩︎D. J. Griffiths. “Quantum Mechanics in Three Dimensions”. In:

*Introduction to Quantum Mechanics*. 2nd. ed. Pearson Printice Hall, 2005. Chap. 4. isbn: 978-0-13-111892-8. ↩︎J. D. Jackson. “Boundary-Value Problems in Electrostatics: II”. In:

*Classical Electrodynamics*. 3rd ed. John Wily & Sons, Inc., 1999. Chap. 3. isbn: 978-0-471-30932-1. ↩︎The mathemeticians' convention typically reverses the role of

`\(\theta\)`

and`\(\phi\)`

so that`\(\theta\)`

measure the azimuth (consistent with polar coordinates in the plane) and`\(\phi\)`

the polar angle, so care must be taken when in interpreting expressions across multiple sources. For more information, see the articles on spherical coordinates at Wolfram MathWorld or Wikipedia. ↩︎