# The Challenge Of Data Analysis For Future CMB Observations

## Introduction

The Cosmic Microwave Background is the most distant observation of photons we can ever make. Last scattered only 300,000 years after the Big Bang it provides a unique picture of the state of the universe at that time. In particular, fluctuations in the CMB directly trace the primordial density perturbations and so provide a powerful discriminant between alternative cosmologies of the very early universe. As a result the search for anisotropies in the CMB has been a cornerstone of cosmology for the last 30 years.

Finally measured by the COBE satellite, the anisotropies proved to be of the order of only one part in a million on a 3K background whose uniformity was otherwise only broken by a dipole induced by the peculiar velocity of the galaxy of the order of one part in a thousand. Despite the tiny scale of these fluctuations, advances in detector technologies have enabled us to consider measuring them to the extraordinary accuracy and resolution necessary to determine the fundamental parameters of cosmology to better than 1% PARAMETERS .

Such measurements include those of the MAXIMA and BOOMERanG projects — described in detail by Lee et al, and Masi et al and de Bernardis et al elsewhere in these Proceedings. These balloon-borne observations have already produced datasets an order of magnitude larger than their predecessors, and in subsequent flights will at least double this size. Beyond this, the MAP and PLANCK satellite missions will yield datasets 1-2 orders of magnitude larger again. The shear size of these datasets makes their analysis a serious computational challenge. It is this challenge, and the current status of our attempts to address it, that are discussed here.

For simplicity we only consider the highly idealised case of extracting a power spectrum of $`𝒩_l`$ multipoles in $`𝒩_b`$ bins from a map of $`𝒩_p`$ pixels obtained from a single time-ordered sequence of $`𝒩_t`$ observations of the sky, the data only comprising CMB signal and Gaussian noise. In practice there are many additional sources of non-Gaussian contamination (in particular both galactic and extra-galactic foreground sources) making observations necessary at a range of frequencies to allow for their subtraction.

## From The Time-Ordered Data To The Map

### Formalism

Our first step is to translate the observation from the temporal to the spatial domain — to make a map MAP  (see also Jaffe et al elsewhere in these Proceedings). Knowing where the detector was pointing at each observation, $`(\theta _t,\psi _t)`$, and adopting a particular pixelization of the sky, we can construct a pointing matrix $`A_{tp}`$ whose entries give the weight of pixel $`p`$ in observation $`t`$. For scanning experiments such as MAXIMA and BOOMERanG this has a particularly simple form

$$A_{tp}=\{\begin{array}{cc}1\hfill & \mathrm{if}(\theta _t,\psi _t)p\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$
(1)   
while other observing strategies would give a more complex structure. The data vector can now be written

$$d=As+n$$
(2)   
in terms of the pixelised CMB signal $`s`$ and time-stream noise $`n`$.

Under the assumption of Gaussianity, the noise probability distribution is

$$P(n)=(2\pi )^{𝒩_t/2}\mathrm{exp}\left\{\frac{1}{2}\left(n^TN^1n+Tr\left[\mathrm{ln}N\right]\right)\right\}$$
(3)   
where $`𝒩`$ is the time-time noise correlation matrix given by

$$Nnn^T$$
(4)   
We can now use equation (2) to substitute for the noise in equation (3), giving the probability of the data for a particular CMB signal as

$$P(d|s)=(2\pi )^{𝒩_t/2}\mathrm{exp}\left\{\frac{1}{2}\left((dAs)^TN^1(dAs)+Tr\left[\mathrm{ln}N\right]\right)\right\}$$
(5)   
Assuming that all CMB maps are a priori equally likely, this is proportional to the likelihood of the signal given the data, and maximizing over $`s`$ gives the maximum likelihood map $`m`$

$$m=\left(A^TN^1A\right)^1A^TN^1d$$
(6)   
Substituting back for the time-ordered data in this map we recover the obvious fact that it is the sum of the true CMB signal and some pixelized noise

$`m`$ $`=`$ $`\left(A^TN^1A\right)^1A^TN^1(As+n)`$ (7)   
$`=`$ $`s+\nu `$   
where this pixel noise

$$\nu =\left(A^TN^1A\right)^1A^TN^1n$$
(8)   
has correlations given by

$`\mathrm{{\rm Y}}`$ $`=`$ $`\nu \nu ^T`$ (9)   
$`=`$ $`\left(A^TN^1A\right)^1`$   
### Computational Requirements

Making the map requires solving equation (6) which is conveniently divided into three steps:

$`\mathrm{{\rm Y}}^1`$ $`=`$ $`A^TN^1A`$   
$`z`$ $`=`$ $`A^TN^1d`$   
$`\mathrm{and}m`$ $`=`$ $`(\mathrm{{\rm Y}}^1)^1z`$ (10)   
The first half of table 1 shows the computational cost of a brute force approach to each of these steps (recall that multiplying an $`[a\times b]`$ matrix and an $`[b\times c]`$ matrix involves $`2\times a\times b\times c`$ operations). Thus for datasets such as MAXIMA-1 or BOOMERanG North America, with $`𝒩_t2\times 10^6`$ and $`𝒩_p3\times 10^4`$, making the map would require of the order of 16 Tb of disc space (storing data in 4-byte precision), 7 Gb RAM<sup>1</sup><sup>1</sup>1Although it is possible to use out-of-core algorithms for operations such as matrix inversion the associated time overhead would be prohibative. We therefore assume that all such operations are carried out in core. (doing all calculations in 8-byte precision) and $`2.4\times 10^{17}`$ operations. If we were able to use a fast 600 MHz CPU at 100% efficiency this would still correspond to over 12 years of run time.

Fortunately there are two crucial structural features to be exploited here. As noted above the pointing matrix $`A`$ is usually very sparse, with only $`𝒩_\alpha `$ non-zero entries in each row. For simple scanning strategies such as MAXIMA, BOOMERanG and PLANCK, $`𝒩_\alpha =1`$, with a single 1 in the column corresponding to the pixel being observed. For a differencing experiment such as COBE or MAP, $`𝒩_\alpha =2`$, with a $`\pm 1`$ pair in the columns corresponding to the pixel pair being observed. Moreover, the inverse time-time noise correlations are (by fiat) both stationary and fall to zero beyond some time-separation much shorter than the duration of the observation

$`N_{tt^{}}^1`$ $`=`$ $`f(|tt^{}|)`$ (11)   
$`=`$ $`0|tt^{}|>\tau 𝒩_t`$   
so that the inverse time-time noise correlation matrix is symmetric and band-diagonal, with bandwith $`𝒩_\tau =2\tau +1`$. The second half of table 1 shows the impact of exploiting this structure on the cost of each step. The limiting step is now no longer constructing the inverse pixel-pixel noise correlation matrix but solving for the map, which is unaffected by these features. For the same datasets making the map now takes of the order of 3.6 Gb of disc, 7 Gb of RAM, and $`7\times 10^{13}`$ operations, or 32 hours of the same CPU time.

Further acceleration of the map-making algorithm must therefore focus on a faster solution the final step, inverting the inverse pixel-pixel noise covariance matrix $`\mathrm{{\rm Y}}^1`$ to obtain the map. However, as we shall see below, even this is not the limiting step overall in current algorithms.

## From The Map To The Power Spectrum

### Formalism

We now want to move to a realm where the CMB observation can be compared with the predictions of various cosmological theories — typically the angular power spectrum. We decompose the CMB signal at each pixel in spherical harmonics

$$s_p=\underset{lm}{}a_{lm}B_lY_{lm}(\theta _p,\psi _p)$$
(12)   
where $`B`$ is the pattern of the observation beam (assumed to be circularly symmetric) in $`l`$-space. The correlations between such signals then become

$$S_{pp^{}}s_ps_p^{}=\underset{lm}{}\underset{l^{}m^{}}{}a_{lm}a_{l^{}m^{}}B_lB_l^{}Y_{lm}(\theta _p,\psi _p)Y_{l^{}m^{}}(\theta _p^{},\psi _p^{})$$
(13)   
For isotropic fluctuations the correlations depend only on the angular separation

$$a_{lm}a_{l^{}m^{}}=C_l\delta _{ll^{}}\delta _{mm^{}}$$
(14)   
and the pixel-pixel signal correlation matrix becomes

$$S_{pp^{}}=\underset{l}{}\frac{2l+1}{4\pi }C_lP_l(\chi _{pp^{}})$$
(15)   
where $`P_l`$ is the Legendre polynomial and $`\chi _{pp^{}}`$ the angle between the pixel pair $`p,p^{}`$. These $`C_l`$ multipole powers completely characterise a Gaussian CMB, and are an otherwise model-independent basis in which to compare theory with observations. In general, due to incomplete sky coverage and low signal-to-noise, we are unable to extract each multipole moment independently. We therefore group the multipoles in bins, adopting a particular spectral shape function $`C_l^s`$ and characterising the CMB signal by its bin powers $`C_b`$ with

$$C_l=C_{b:lb}C_l^s$$
(16)   
Since the signal and noise are assumed to be realisations of independent Gaussian processes the pixel-pixel map correlations are

$`M_{pp^{}}`$ $``$ $`mm_T`$ (17)   
$`=`$ $`ss_T+\nu \nu ^T`$   
$`=`$ $`S+\mathrm{{\rm Y}}`$   
and the probability distribution of the map given a particular power spectrum $`𝒞`$ is now

$$P(m|𝒞)=(2\pi )^{𝒩_p/2}\mathrm{exp}\left\{\frac{1}{2}\left(m^TM^1m+Tr\left[\mathrm{ln}M\right]\right)\right\}$$
(18)   
Assuming a uniform prioir for the spectra, this is proportional to the likelihood of the power spectrum given the map. Maximizing this over $`𝒞`$ then gives us the required result, namely the most likely CMB power spectrum underlying the original observation $`d`$.

Finding the maximum of the likelihood function of equation (18) is generically a much harder problem than making the map. Since there is no closed-form solution corresponding to equation (6) we must find both a fast way to evaluate the likelihood function at a point, and an efficient way to search the $`𝒩_b`$-dimensional parameter space for the peak. The fastest general method extant is to use Newton-Raphson iteration to find the zero of the derivative of the logarithm of the likelihood function MLA . If the log likelihood function

$$(𝒞)=\frac{1}{2}\left(m^TM^1m+Tr\left[\mathrm{ln}M\right]\right)$$
(19)   
were quadratic, then starting from some initial guess at the maximum likelihood power spectrum $`𝒞_o`$ the correction $`\delta 𝒞_o`$ that would take us to the true peak would simply be

$$\delta 𝒞_o=\left(\left[\frac{^2}{𝒞^2}\right]^1\frac{}{𝒞}\right)_{𝒞=𝒞_o}$$
(20)   
Since the log likelihood function is not quadratic, we now take

$$𝒞_1=𝒞_o+\delta 𝒞_o$$
(21)   
and iterate until $`\delta 𝒞_n0`$ to the desired accuracy. Because any function is approximately quadratic near a peak, if we start searching sufficiently close to a peak this algorithm will converge to it. Of course there is no guarantee that it will be the global maximum, and in general there is no certainty about what ‘sufficiently close’ means in practice. However experience to date suggests that the log likelihood function is sufficiently strongly singley peaked to allow us to use this algorithm with some confidence.

### Computational Requirements

The core of the algorithm is then to calculate the first two derivatives of the log likelihood function with respect to the multipole bin powers

$`{\displaystyle \frac{}{𝒞_b}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m^TM^1{\displaystyle \frac{S}{𝒞_b}}M^1mTr\left[M^1{\displaystyle \frac{S}{𝒞_b}}\right]\right)`$ (22)   
$`{\displaystyle \frac{^2}{𝒞_b𝒞_b^{}}}`$ $`=`$ $`m^TM^1{\displaystyle \frac{S}{𝒞_b}}M^1{\displaystyle \frac{S}{𝒞_b^{}}}M^1m+{\displaystyle \frac{1}{2}}Tr\left[M^1{\displaystyle \frac{S}{𝒞_b}}M^1{\displaystyle \frac{S}{𝒞_b^{}}}\right]`$ (23)   
Evaluating these derivatives comes down to solving the $`𝒩_b𝒩_p+1`$ linear systems

$`z`$ $`=`$ $`M^1m`$ (24)   
$`\mathrm{and}W_b`$ $`=`$ $`M^1{\displaystyle \frac{S}{𝒞_b}}b`$ (25)   
and assembling the results

$`{\displaystyle \frac{}{𝒞_b}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m^TW_bzTr\left[W_b\right]\right)`$   
$`{\displaystyle \frac{^2}{𝒞_b𝒞_b^{}}}`$ $`=`$ $`m^TW_bW_b^{}z+{\displaystyle \frac{1}{2}}Tr\left[W_bW_b^{}\right]`$ (26)   
Table 2 shows the computational cost of these operations, where solving the linear systems has been optimised by first Cholesky decomposing the matrix $`M`$. Solving equation (20) has been excluded since its cost is entirely negligible, depending only on the number of multipole bins $`𝒩_b𝒩_p`$. Obtaining the maximum likelihood power spectrum for the same datasets as above, with $`𝒩_p3\times 10^4`$ and $`𝒩_b20`$, then requires of the order of 150 Gb disc, 14 Gb of RAM, and $`10^{15}`$ operations per iteration, or 21 days of our 600 MHz CPU time.

Such numbers are at least conceivable; however, as shown in table 3, the scaling with map size pushes forthcoming balloon observations well beyond the capacity of the most powerful single processor machine — and even allowing for the continued doubling of computer power every 18 months predicted by Moore’s ‘law’ we would still have to wait 20 years for a serial machine capable of handling the BOOMERanG Long Duration Balloon flight data. Moreover, these timings are for a single iteration (and typically the alogorithm needs $`𝒪(10)`$ iterations to converge) for a single run through the dataset, although undoubtedly we will want to perform several slightly different runs to check the robustness of our analysis.

One way to increase our capability now is to move to parallel machines, such as the 640-processor Cray T3E at NERSC. Since the algorithm is limited by matrix-matrix operations (Cholesky decomposition and triangular solves) we are able to exploit the most optimised protocols — the level 3 BLAS — and the DEC Alpha chips’ capacity to perform an add and a multiply in a single clock cycle. Coupled with a finely-tuned dense packing of the matrices on the processors this has enabled us to sustain upwards of 2/3 of the theoretical peak performance of 900MHz. This enables us to increase the limiting datasize to around 80,000 pixels.

## Future Prospects

We have seen that existing algorithms are capable of dealing with CMB datasets with at most $`10^5`$ pixels. Over the next 10 years a range of observations are expected to produce datasets of $`5\times 10^5`$ (BOOMERanG LDB), $`10^6`$ (MAP) and $`10^7`$ (PLANCK) pixels that will necessarily require new techniques. This is an ongoing area of research in which some progress has been made in particular special cases.

The limiting steps in the above analysis are associated with operations involving the pixel-pixel correlation matrices for the noise $`\mathrm{{\rm Y}}`$, the signal $`S`$, and most particularly their sum $`M`$. The problem here is the noise and the signal have different natural bases. The inverse noise correlations are symmetric, band-diagonal and approximately circulant in the time domain, while the signal correlations are diagonal in the spherical harmonic domain. However there is no known basis in which the map correlations take a similarly simple form.

If the noise is uncorrelated between pixels and is azimuthally symmetric — as has been argued will be the case for the MAP satellite due to its chopped observing strategy — then the pixel-pixel data correlation matrix can be dramatically simplified, reducing the analysis to $`𝒪(𝒩_p^{3/2})`$ in storage and $`𝒪(𝒩_p^2)`$ in operations APPROX . Although some work has also been done extending this to observations with marginal azimuthal asymmetry it is still inapplicable for the spatially correlated noise inherent to the simple scanning strategies adopted by MAXIMA and BOOMERanG (which also face the problem of irregular sky coverage) or PLANCK; for such observations the problem remains unsolved.

## Acknowledgements

The author would like to thank the organisers of the Rome 3K Cosmology Euroconference for the opportunity to participate in an excellent meeting, and Andrew Jaffe, Pedro Ferriera, Shaul Hanany, Xiaoye Li, Osni Marques, and Radek Stompor for many useful discussions.