MADmap: A Fast Parallel Maximum Likelihood CMB Map Making Code

Introduction

MADmap is a parallel code which will produce maximum likelihood sky maps from time ordered instrument data. This code was written with large data sets in mind, but works well on any time ordered data set which has piecewise stationary noise. MADmap is part of MADCAP++ which also includes a power spectrum estimator (the full version of MADCAP++ has not yet been released). MADCAP++ is a version of MADCAP intended for large data sets. For more information about MADCAP see Julian Borrill's web page .

Data Model

This code is built around the following data model:
where d is the time ordered data set n is a time domain vector of piece wise stationary Gaussian noise, s is the signal in some basis other than the time domain, and A is the pointing matrix which is the linear operator which projects from the signal basis to the time domain. This is a very general data model. Many different experimental data sets can be formulated this way. Once the data is formulated in this way the maximum likelihood map can be derived with MADmap.

Maximum Likelihood Map Equation

MADmap solves the map making equation:
A is the n_t by n_p pointing matrix where n_t is the number of time samples and n_p is the number of pixels. N is the time time noise covariance matrix. d is the time ordered data set (the observation with noise). m is the pixel domain maximum likelihood estimate of the noiseless signal map given N and d.

Method of Calculation

This calculation is done by assuming that inv(N) is band diagonal and piecewise Toeplitz (note that because N is a covariance matrix, N and inv(N) are symmetric and positive definite). That is to say that we assume that the operation of inv(N) is equivalent to a set of convolutions with band limited kernels (i.e. kernels which are zero beyond some band width that is shorter than the length of the observation). This convolution is done in the Fourier domain by way of the convolution theorem, so the operation of inv(N) can be done in O(n_t ln( n_b )) calculations where n_b is the band width of inv(N)'s diagonal (i.e. the width of the kernels).

Given that we can calculate the operation of inv(N) on a vector, we still need to compute the operation:

This is done by preconditioned conjugate gradient method, where each iteration of the conjugate gradient routine involves the inv(N) operator which is calculated as described above. The preconditioner is the pixel domain diagonal matrix with the one over the number of observations in each pixel along the diagonal.

MADmap on seaborg.nersc.gov

MADmap is installed and publicly available on the NERSC IBM-SP called seaborg. The directory is:

~cmc/MADmap

and the executable is

~cmc/MADmap/MADmap