MADmap evaluates the following equation:
m = inv(A'*inv(N)*A)*A'*inv(N)*d
where A is the pointing matrix (and A' is A transpose), N is the time time noise covariance matrix (and inv(N) is the inverse of N), d is the time ordered data vector, and m is the maximum likelihood pixel map vector given N and d.
This calculation is done by assuming that inv(N) is band diagonal and 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 convolution with a band limited kernel (i.e. a kernel which is zero beyond some band width which 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 kernel).
Given that we can calculate the operation of inv(N) on a vector, we still need to evaluate the following linear operator:
inv(A'*inv(N)*A)
This is done by preconditioned conjugate gradient method, where each iteration of the conjugate gradient routine involves the inv(N) operator which we calculate 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. For most CMB data sets this converges with a relative residual of 1e-5 in approximately 30 iterations.
The command line arguments are as follows:
| (A'*inv(N)) * (A*m - d) |
e = ---------------------------
| (A'*inv(N)) * d |
Note that * is the dot product, x' means x transpose, and | x | is the Euclidean norm of the vector x ( i.e. sqrt(x'*x) ).
Given these definitions MADmap should scale as
O( i*n*( f*ln(f)/(f-b) + m ) / p )
in time, and
O( n )
There are two file types used by this code. One file type is used to store the time ordered data, the other is used to store the inverse time time noise correlations. Both of these files types were derived from the original MADCAP file types, but have been modified slightly. The format of the two file types are described below:
Time ordered data file
This is a binary file which contains the time ordered data. This file records the data collected by an instrument and the pointing information in a time series. Such a file starts with a header composed of four eight byte integers. In order, these are (a) the first time index, (b) the last time index, (c) the number of pixels observed during each time sample, and finally the total number of pixels. For each time sample there is set of data recorded in the file. If there are n pixels observed during each time sample, then there are 2*n+1 numbers recorded for each time sample. The first number is the observation stored as an eight byte floating point number. This is followed by n pairs of numbers. These pairs are the weight of the observation and the pixel index this weight applies to. The weights are stored as four byte floating point numbers, and the pixel indexes are stored as four byte integers. If the pixel number is negetive then the weight is effectively zero (i.e. fewer pixels were observed during time samples which have negative pixel indexes).
In summary the body of the time ordered data file could
be described as follows:
for each observation
data (8-byte float)
for each pixel observed
weight (4-byte float)
pixel index (4-byte int)
Inverse time-time noise correlation files
These are binary files which describe the inverse time time noise correlation matrix. There is a file for each stationary time interval of the data. These files should start with the same prefix and end with the suffixes ".0", ".1", ".2", etc.
MADmap makes the assumption that the inverse of the time time noise correlation matrix is block diagonal, with one block for each stationary time interval in the data. These blocks are assumed to be Toeplitz and band diagonal (the blocks are also symmetric because this is a correlation matrix). This is not generally true, but with care the noise can be adequately characterized this way. If we know the non-zero elements of the first row of each Toeplitz block, and where along the diagonal each block resides, we can completely describe the inverse time time noise correlation matrix. All of this information is given by the inverse time time noise correlation files.
One file describes each block of the matrix. These files start with a header composed of three eight byte integers. In order, these are (a) the number of non-zero elements in the first row of this block minus one (i.e. the number of correlations), (b) the time index of the upper left hand corner of the block, (c) the time index of the lower right hand corner of the block. This header is followed by a series of eight byte floating point numbers. These are the non-zero elements of the first row of the block which covers the time samples given in the header. Note that if the number of correlations given by the header were c, then there would be c+1 floating point numbers in the file and this would correspond to a band width of 2*c+1. When running the code a band width is specified; this is a maximum band width allowed. If the number of correlations in a file corresponds to a larger band width than the one specified at the command line, then the extra correlations in the file are ignored, and the computation is done assuming that they are zero.
Any region of the time stream which is not within a block described by one of the inverse time time noise correlation files will be ignored (i.e. these sections of the inverse noise covariance matrix are assumed to be zero).
The MLE map
The primary output of the executable is the maximum likelihood map. The name of this file is specified by the user. As mentioned in the description, this is an eight byte floating point raw binary file.
Diagnostic files
Two diagnostic files will be created by this code. The first is a file which contains the number of observations in each pixel (as eight byte floating point numbers). This is actually just the pointing matrix doted with a vector of ones:
A'*1
This file will have the same name as the time ordered data file, but it will end in ".obsDensity" (i.e. if the time ordered data file were called "t_data" then the created file would be called "t_data.obsDensity". The other file to be created will be called "mapFirstStep.dat" This is the result of the first computation:
A'*inv(N)*d.
This is the so called "noise weighted map" which is stored as an eight byte floating point binary file.
More extensive html documentation should have been distributed with this software. This html documentation can also be found on the web at http://www.nersc.gov/~cmc/MADmap/doc
MADmap is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
Christopher Cantalupo <cmc@nersc.gov>
Send bug reports or comments to the above address.