text:
# MAPCUMBA : a fast iterative multi-grid map-making algorithm for CMB experiments

## 1 Introduction

As cosmology enters the era of “precision”, it enters simultaneously the era of massive data sets. This has in turn showed the need for new data processing algorithm. Present and future CMB experiments in particular face some new and interesting challenges Bond et al. (1999). If we accept the now classical point of view of a four steps data analysis pipeline : i/ from *time-ordered data* (TOD) to maps of the sky at a given frequency, ii/ from frequency maps to (among others) a temperature map, iii) from a temperature map to its power spectrum $`C_{\mathrm{}}`$, iv/ from the power spectrum to cosmological parameters and characteristics of the primordial spectrum of fluctuation, the ultimate quantities to be measured in a given model. The work we are presenting focus on the first of these issues, namely the map-making step.

Up to the most recent experiments, maps could be made by a brute force approach amounting to solve directly a large linear problem by direct matrix inversion. Nevertheless the size of the problem, and the required computing power, grows as a power law of the data set size, and the limits of this kind of method have now been reached Borrill (2000). Whereas the most efficient development in this massive computing approach, i.e. the MADCAP package Borrill (1999) has been applied to the recent BOOMERanG and MAXIMA experiments de Bernardis et al. (2000); Hanany et al. (2000) some faster and less consuming solutions based on iterative inversion algorithms have now been developed and applied too Prunet et al. (2000). This is in the same spirit as Wright et al. (1996), and we definitely follow this latter approach.

Design goals are to use memory sparingly by handling only columns or rows instead of full matrixes and to increase speed by minimising the number of iterations required to reach the sought convergence of the iterative scheme. We fulfill these goals by an iterative multi-grid Jacobi algorithm. As recalled below, an optimal method involves using the noise power spectrum. We have thus investigated the possibility of a joint noise (and signal) evaluation using this algorithm.

Section 2 presents in detail the method and its implementation, while section 3 demonstrates its capabilities by using simulations of two on-going balloon CMB experiments, ARCHEOPS111http://www-crtbt.polycnrs-gre.fr/archeops/Benoit et al. (2000) and TopHat222http://topweb.gsfc.nasa.gov/. The results are discussed in section 4, together with the problem of the evaluation of the noise properties, as well as possible extensions.

## 2 The method

### 2.1 Optimal map making

The relation between the sky map we seek and the observed data stream may be cast as a linear algebra system Wright et al. (1996); Tegmark (1997). Let t and p indices denote quantities in the temporal and spatial domains, and group as a data vector, $`d_t`$, and a noise vector the temporal stream of collected data and the detector noise stream, both of dimension $`𝒩_{tod}`$. We then have

$$d_t=A_{tp}x_p+n_t,$$
(1)
where $`A_{tp}x_p`$ is the signal vector given by the observation of the unknown pixelised sky map, $`x_p`$, which has been arranged as a vector of dimension $`𝒩_{pix}`$. The $`𝒩_{tod}\times 𝒩_{pix}`$ “observation” matrix $`A`$ therefore encompasses the scanning strategy and the beam pattern of the detector.

In the following, we restrict to the case when the beam pattern is symmetrical. We can therefore take $`x_p`$ to be a map of the sky which has already been convolved with the beam pattern, and $`A`$ only describes how this sky is being scanned. For the total power measurement (i.e. non-differential measurement) we are interested in here, the observation matrix $`A`$ then has a single non-zero element per row, which can be set to one if $`d`$ and $`x`$ are expressed in the same units. The model of the measurement is then quite transparent: each temporal datum is the sum of the pixel value to which the detector is pointing plus the detector noise.

The map-making step then amounts to best solve for $`x`$ given $`d`$ (and some properties of the noise). We shall seek a *linear* estimator of $`x_p`$,

$$\widehat{x}_p=W_{pt}d_t.$$
(2)
To motivate a particular choice of the $`𝒩_{pix}\times 𝒩_{tod}`$ matrix $`W`$, a Bayesian approach is convenient. Indeed we are seeking the optimal solution to this inversion problem which maximises the probability of a deduced set of theory parameters (here the map $`x_p`$) given our data ($`d_t`$) by maximising $`𝒫(x|d)`$. Bayes’ theorem simply states that

$$𝒫(x|d)=\frac{𝒫(d|x)𝒫(x)}{𝒫(d)}.$$
(3)
If we *do not assume any theoretical prior*, then $`x`$ follows a uniform distribution as well as $`d`$. Therefore,

$`𝒫(x|d)`$ $``$ $`𝒫(d|x).`$ (4)
If we further assume that *the noise is Gaussian*, we can write

$`𝒫(x|d)`$ $``$ $`exp(n_t^TN_{tt^{}}^1n_t^{}/2)`$ (5)
$``$ $`exp((dAx)_t^TN_{tt^{}}^1(dAx)_t^{}/2)`$ (6)
$``$ $`exp(\chi ^2/2)`$ (7)
where $`N_{tt^{}}=<nn^T>_{tt^{}}`$ is the noise covariance matrix. In this particular case, maximising $`𝒫(x|y)`$ amounts to find the least square solution which was used to analyse the “COBE” data Jansen and Gulkis (1992),

$$W=\left[A^TN^1A\right]^1A^TN^1.$$
(8)
In this paper we will actually deal only with this estimator. Nevertheless as a next iteration in the analysis process, we could incorporate various theoretical priors by expliciting $`𝒫(x)`$. For example, it is often assumed *a Gaussian prior* for the theory, i.e. $`𝒫(x)exp(x_p^TC_{pp^{}}^1x_p^{}/2)`$ where $`C_{pp^{}}=x_px_p^{}^T`$ is the signal covariance matrix. In that case the particular solution turns out to be the Wiener filtering solution Zaroubi et al. (1995); Bouchet and Gispert (1996); Tegmark and Efstathiou (1996); Bouchet and Gispert (1998):

$$W=\left[C^1+A^TN^1A\right]^1A^TN^1.$$
(9)
But this solution may always be obtained by a further (Wiener) filtering of the COBE solution, and we do not consider it further.

The prior-less solution demonstrates that as long as the (Gaussian) instrumental noise is not white, a simple averaging (co-addition) of all the data points corresponding to a given sky pixel is not optimal. If the noise exhibits some temporal correlations, as induced for instance by a low-frequency $`1/f`$ behavior of the noise spectrum which prevails in most CMB experiments, one has to take into account the full time correlation structure of the noise. Additionally this expression demonstrates that even if the noise has a simple time structure, the scanning strategy generically induces a non-trivial correlation matrix $`\left[A^TN^1A\right]^1`$ of the noise map.

Even if the problem is well posed formally, a quick look at the orders of magnitude shows that the actual finding of a solution is non trivial task. Indeed a brute force method aiming at inverting the full matrix $`\left[A^TN^1A\right]^1`$, an operation scaling as $`𝒪(𝒩_{pix}^3)`$, is already hardly tractable for present long duration balloon flights as MAXIMA, BOOMERanG, ARCHEOPS or TopHat where $`𝒩_{tod}10^6`$ and $`𝒩_{pix}10^5`$. It appears totally impractical for PLANCK since for a single detector (amid $`10`$s) $`𝒩_{tod}10^9`$ and $`𝒩_{pix}10^7`$!

One possibility may be to take advantage of specific scanning strategies, and actually solve the inverse of the convolution problem as detailed in Wandelt and Gorski (2000). This amounts to deduce the map coefficients $`a_{lm}`$ in the spherical harmonic basis through a rotation of a Fourier decomposition of the observed data. The map will then be a simple visualisation device, while the $`a_{lm}`$ would be ready to use directly for a component separation (as in Bouchet and Gispert (1996); Tegmark and Efstathiou (1996); Bouchet and Gispert (1998)) and the CMB power spectrum estimate. While potentially very interesting, this approach will not be generally applicable (at least efficiently), and we now turn to a practical (general) solution of equation (8) by iterative means.

### 2.2 Practical implementation for large data sets

We solve the map-making problem by adapting to our particular case the general “multi-grid method” Press et al. (1992). Multi-grid methods are commonly used to speed up the convergence of a traditional relaxation method (in our case the Jacobi method, as in Prunet et al. (2000)) defined at resolution $`\mathrm{}_{max}`$ (see below). A set of recursively defined coarser grids ($`\mathrm{}<\mathrm{}_{max}`$) are used as temporary computational space, in order to increase the convergence rate of the relaxation process. To be fully profitable, this algorithm implies for each resolution both a rebinning in space (resolution change) and in time (resampling).

In this paper, we use the HEALPix pixelisation of the sphere Górski et al. (1998). In this scheme, the sphere is covered by 12 basic quadrilaterals, further divided recursively into pixels of equal area. The map resolution is labeled by $`N_{side}`$: the number of pixels along the side of one basic quadrilateral. Hence, $`N_{side}=1`$ means that the sphere is covered by 12 large pixels only. The number of pixels is given by $`N_{pix}=12N_{side}^2`$. $`N_{side}=256`$ corresponds to a pixel size of $`13.7`$ arcmin. For practical reasons, we need to define the resolution $`k`$ of a HEALPix map as

$$N_{side}=2^k$$
(10)
The “nested”pixel numbering scheme of HEALPix Górski et al. (1998) allows an easy implementation of the coarsening ($`kk1`$) and refining ($`kk+1`$) operators that we use intensively in our multigrid scheme.

Let us now get into the details of our implementation and discuss successively the exact system we solve, the way we solve it and the actual steps of the multi-grid algorithm.

#### 2.2.1 Determining the working resolution $`k_{max}`$

We aim at solving for the optimal map $`\widehat{x}_k`$ at a given spatial resolution $`k`$ using

$$A_k^TN^1A_k\widehat{x}_k=A_k^TN^1d,$$
(11)
where $`A_k`$ is the “observation” operator (from spatial to temporal domain) and $`A_k^T`$ is the “projection” operator (from temporal to spatial domain). In a noise-free experiment, the optimal map would be straightforwardly given by the co-added map (introducing the “co-addition” operator $`P_k`$)

$$\widehat{x}_k=P_kd(A_k^TA_k)^1A_k^Td$$
(12)
The time line is given by $`d=A_{\mathrm{}}x_{\mathrm{}}+n`$ where $`x_{\mathrm{}}`$ is the sky map at “infinite” resolution ($`k=+\mathrm{}`$ in our notations). In order to check the accuracy of this trivial noise-free map making, it is natural to compute the residual in the time domain with $`n=0`$

$$p_k=A_k\widehat{x}_kd=A_k\widehat{x}_kA_{\mathrm{}}x_{\mathrm{}}$$
(13)
which will be non-zero in practice, as soon as one works with finite spatial resolution. We call this residual the pixelisation noise. Since we assume here that the instrumental beam is symmetric, the sky map is considered as the true sky convolved by, say, a Gaussian beam of angular diameter $`\mathrm{\Delta }\theta _B`$. This introduces a low-pass spatial filter in the problem. In other words, as the resolution increase, the pixelisation noise should decrease towards zero. We have estimated that the order of magnitude of the pixelisation noise can be approximated by

$$p_kx_{\mathrm{}}\frac{\mathrm{\Delta }\theta _k}{\mathrm{\Delta }\theta _B}$$
(14)
The norms used in the above formula can be either the maximum over the time line (a very strong constraint) or the variance over the time line (a weaker constraint). Since the pixelisation noise is strongly correlated with the sky signal, point sources or Galaxy crossings are potential candidates for large and localised bursts of pixelisation noise. The correct working resolution $`k_{max}`$ is set by requiring that the pixelisation noise remains low compared to the actual instrumental noise, or equivalently

$$\mathrm{\Delta }\theta _{k_{max}}\frac{\mathrm{\Delta }\theta _B}{S/N}$$
(15)
Most of the CMB experiments are noise dominated along the time line, constraining the effective map resolution to be of the order of the instrumental beam or even larger. Note however that the pixelisation noise is strongly non Gaussian (point sources or Galaxy crossings) and can be always considered as a potential source of residual stripes in the final maps.

#### 2.2.2 The basic relaxation method: an approximate Jacobi solver

Instead of solving for $`\widehat{x}`$ we perform the change of variable

$$\widehat{y}=\widehat{x}Pd$$
(16)
and solve instead for $`\widehat{y}`$ which obeys

$$PN^1A\widehat{y}=PN^1(dAPd)\text{or}M\widehat{y}=b$$
(17)
where we have multiplied each side of equation (11) by $`(A^TA)^1`$, the pixel hit counter. From now on, we also assume that the noise in the timestream is stationary and that its covariance matrix is normalised so that $`\text{diag}N_{tt^{}}^1=I`$. The previous change of variable allows us to subtract the sky signal from the data: since we have chosen a resolution high enough to neglect the pixelisation noise, we have indeed $`APdA_{\mathrm{}}x_{\mathrm{}}+APn`$ and, consequently, $`dAPdnAPn`$. The map making consists in two step: first compute a simple co-added map from the time line, and second, solve equation (17) for the stripes map $`\widehat{y}`$. The final optimal map is obtained by adding these two maps.

It is worth mentioning that the stripes map is completely independent of the sky signal, as soon as the pixelisation noise can be ignored. It depends only on the scanning strategy through the matrix $`A`$ and on the noise characteristics through the matrix $`N`$. Even if in principle this change of variable is irrelevant since it does not change the matrix to be inverted, it does in practice since $`dAPd`$ is free from the fast variations of $`d`$ (up to the pixelisation noise), as e.g. the Galaxy crossings, which are numerically damaging to the quality of the inversion.

To solve for equation 17, we follow the ideas of Prunet et al. (2000) and apply an approximate Jacobi relaxation scheme. The Jacobi method consists in the following relaxation step to solve for our problem

$$\widehat{y}^{n+1}=R\widehat{y}^n+D^1b\text{and}\widehat{y}^0=0$$
(18)
where $`D`$ is the diagonal part of the matrix $`M`$ and $`R`$ is the residual matrix $`R=ID^1M`$. Computing the diagonal elements of $`M=PN^1A`$ is rather prohibitive. The idea of Prunet et al. (2000) is to approximate $`DI`$ by neglecting the off-diagonal elements of $`N^1`$. The residual matrix then simplifies greatly and writes

$$R=P(IN^1)A$$
(19)
The approximate Jacobi relaxation step is therefore defined as

$$\widehat{y}^{n+1}=R\widehat{y}^n+b\text{and}\widehat{y}^0=0$$
(20)
One clearly sees that if this iterative scheme converges, it is towards the full solution of equation 17. To perform these successive steps, it is extremely fruitful to remember the assumed stationarity of the noise. Indeed whereas this assumption implies a circulant noise covariance matrix in real space, it translates in Fourier space in the diagonality of the noise covariance matrix. This is naturally another formulation of the convolution theorem, since a stationary matrix acts on a vector as a convolution, and a diagonal matrix acts as a simple vector product, thus a convolution in real space is translated as a product in Fourier space. The point is that the manipulation of the matrix $`N^1`$ is considerably lighter and will be henceforth performed in Fourier space. Applying the matrix $`R`$ to a map reduces then to the following operations in order

1. “observe” the previous stripes map $`\widehat{y}^n`$
2. Fourier transform the resulting data stream
3. apply the low-pass filter $`W=IN^1`$
4. inverse Fourier transform back to real space
5. co-add the resulting data stream into the map $`\widehat{y}^{n+1}`$.

Assuming that the normalised noise power spectrum can be approximated by $`P(f)=1+\left(f_0/f\right)^\alpha `$, the low-pass filter associated to each relaxation step is given by

$$W(f)=\frac{f_0^\alpha }{f_0^\alpha +f^\alpha }$$
(21)
Since both $`A`$ and $`P`$ are norm-preserving operators, the norm of the increment $`\mathrm{\Delta }\widehat{y}^n=\widehat{y}^n\widehat{y}^{n1}`$ between step $`n`$ and $`n+1`$ decreases as $`\mathrm{\Delta }\widehat{y}^{n+1}W(f_{min})\mathrm{\Delta }\widehat{y}^n`$, where $`f_{min}`$ is a minimal frequency in the problem. Since $`W(f_{min})<1`$, we see that the approximate Jacobi relaxation scheme will converge in every case, which is good news. On the other hand, since $`W(f_{min})1`$, the actual convergence rate of the scheme is likely to be very, very slow, which is bad news (cf. figure 6 for a graphical illustration) . The fact that this algorithm is robust, but dramatically slow is a well-known property of the Jacobi method. The multi-grid method is also well known to solve this convergence speed problem. Note that if the convergence is reached, the solution we get is the optimal solution, i.e. similar to the one that would be obtained by a full matrix inversion.

#### 2.2.3 Multigrid relaxation

The multi-grid method for solving linear problems is described in great details in Press et al. (1992). At each relaxation step at level $`k=k_{max}`$, our target resolution, we can define the error $`e_k^n=\widehat{y}_k^n\widehat{y}_k`$ and the residual $`r_k^n=M_k\widehat{y}_k^nb_k`$. Both are related through

$$M_ke_k^n=r_k^n$$
(22)
If we are able to solve exactly for equation (22), the overall problem is solved. The idea of the multi-grid algorithm is to solve approximately for equation (22) using a coarser “grid” at resolution $`k1`$, where the relaxation scheme should converge faster. We thus need to define a fine-to-coarse operator, in order to define the new problem on the coarse grid and solve for it. We also need a coarse-to-fine operator in order to inject the solution onto the fine grid. The approximation to the error $`e_k^n`$ is finally added to the solution. The coarse grid solver is usually applied after a few iterations on the fine level have been performed (in practice we perform 3 to 5 iterations). Naturally, since the solution to the problem on the coarse level relies also on the same relaxation scheme, it can be itself accelerated using an even coarser grid. This naturally leads to a recursive approach of the problem.

We defined our fine-to-coarse operator to be an averaging of the values of the 4 fine pixels contained in each coarse pixel. The coarse-to-fine operator is just a straight injection of the value of the coarse pixel into the 4 fine pixels. The most important issue is the temporal rebinning of the data stream since the speed of the iterative scheme at a given level is actually set by the two Fourier transforms. We performed that resampling by simply taking each time we go up one level one point out of two. At the lower resolutions, this reduction is such that the iteration cost is negligible when compared to that at higher $`k`$; it allows fast enough iterations that full convergence may be reached. In practice we choose a minimal level $`k=3`$ defined by $`N_{side}=8`$ and iterate a few hundred times to reach exact convergence.

Finally the navigation through levels allows several options to be taken. Either we go up and down through all the levels successively (the so-called “V-cycle”) or we follow more intricate paths (e.g. the “W-cycle” where at a given level we go down and up all the lower levels before going up and vice-versa). Since it turns out that levels are relatively disconnected in the sense that the scales well solved at a given resolution are only slightly affected by the solution at a higher level, the “V-cycle” is appropriate and is the solution we adopt in the two following configurations.

## 3 Practical application to ARCHEOPS and TopHat experiments

We now aim at demonstrating the capabilities of this algorithm with simulated data whose characteristics are inspired by those of the ARCHEOPS and TopHat experiments.

### 3.1 Simulation

The ARCHEOPS and TopHat experiments are very similar with respect to their scanning strategy. Indeed both use a telescope that simply spins at a constant rate (respectively $`3`$ and $`4`$ RPM) about its vertical axis. Thus due to Earth rotation the sky coverage of both is performed through scan circles whose axis slowly drifts on the sky. Nevertheless, because of the different launch points (respectively on the Arctic circle (Kiruna, Sweden) and in Antarctica (McMurdo)) and their different telescope axis elevation (respectively $`41^o`$ or $`12^o`$) they do not have the same sky coverage.

Otherwise the two experiments do not use the same bolometers technology, neither the same bands nor have the same number of bolometers. But even if we try to be fairly realistic, our goal though is not to compare their respective performances but rather to demonstrate two applications of our technique in different settings. We then simulate for each a data stream of $`24\mathrm{hrs}`$ duration with respectively a sampling frequency of $`171`$ and $`64\mathrm{Hz}`$. The TODs contain realistic CMB and Galactic signal for a band at $`143\mathrm{GHz}`$. Note that this is a one day data stream for TopHat (out of 10 expected) and that this frequency choice is more appropriate for ARCHEOPS than for TopHat (whose equivalent band is around $`156\mathrm{GHz}`$), but this is realistic enough for our purpose. We generated a Gaussian noise time stream with the following power spectrum $`P(f)(1+(f_{knee}/f)^\alpha )^1`$ with $`f_{knee}=0.24`$ and $`1\mathrm{Hz}`$, and $`\alpha =1.68`$ and $`1`$. The noise amplitude per channel is choosen so that it corresponds for ARCHEOPS (24 hours) and TopHat333Lloyd Knox private communication (10 days of flight) to $`30/8\mu \text{K}`$ on average per $`20^{}`$ pixel, with a beam FWHM of $`10^{}`$ / $`20^{}`$.

We introduced $`5`$ distinct levels of resolution defined by their $`N_{side}`$ parameter in the HEALPix package Górski et al. (1998). The higher resolution level is imposed by the pixelisation noise level requirement (section 2.2.3) to $`N_{side}=256`$ (pixel size $`13.7^{}`$) whereas the lower one is $`N_{side}=8`$ (pixel size $`7.3^o`$). Therefore these two configurations each offer an interesting test since they differ by the sky coverage and the noise power spectrum. We iterate $`3`$ times at each level except at the lowest one where we iterate $`100`$ times.

### 3.2 Results

The algorithm is as efficient in both situations. Whereas for ARCHEOPS whose timeline is longer due to the higher sampling frequency it took 2.25 hours on a SGI ORIGIN 2000 single processor, it took around 1.5 hours for the TopHat daily data stream.

In figures 1 and 2 we depict from top to bottom and from left to right the initial co-added data, the reconstructed noise map, the hit map, i.e. the number of hit per pixels at the highest worked out resolution, the initial signal map as well as the reconstructed signal map and the error map. We see that the destriping is excellent in both situations and the signal maps recovered only contain an apparently isotropic noise. We note the natural correlation between the error map and the hit map. Finally, we stress that no previous filtering at all has been applied.

Figure 3 shows how the noise map is reconstructed at various scales. This is a mere illustration of our multi-grid work.

### 3.3 Tests

At this point, some tests are required to cautiously validate our method. First, as was stated below, as soon as the iterative algorithm has converged, the solution is by construction the optimal solution, similar to the one that would be obtained by the full matrix inversion. As a criterium for convergence we required the $`2`$-norm of residuals to be of the order of the machine accuracy.

We initially have a Gaussian random noise stream fully characterised by its power spectrum. Therefore an important test is to check wether the deduced noise stream (by “observing” the stripe map, see § 4.3 for further details definition) is Gaussian and has the same power spectrum. On figure 4 we ensure that the evaluated noise time stream is indeed gaussian. As depicted on figure 5, where we plot in the ARCHEOPS case both the analytical expression of the spectrum according to which we generate the timeline and its recovered form, the agreement is excellent. We recall that we assumed at the beginning a perfect knowledge of the noise in order to define the filters. This is naturally unrealistic but the issue of noise evaluation is discussed in section 4.3 below. We plotted as well the probability distribution function (PDF) of the final error map, i.e. the recovered noisy signal map minus the input signal map (figure 4). This PDF is well fitted by a gaussian whose parameters are given in figure 4. The PDF of the error map displays some non-gaussian wings. Let us recall here that this is no surprise here because of the non-uniform sky coverage as well as the slight residual stripping, both due to the non-ideal scanning strategy, i.e. that produces a non-uniform white noise in pixel space.

Another particularly important test consists in checking the absence of correlation between the recovered noise map and the initial signal map. We could not find any which is no surprise since we are iterating on a noise map (see § 2.2.2) which does not contain any signal up to the pixelisation noise, that is ensured to be negligible with regards to the instrumental noise by choice of the resolution $`k_{max}`$ (see § 2.2.1).

## 4 Discussion

### 4.1 Why is such an algorithm so efficient ?

The efficiency of such a method can be intuitively understood. Indeed, although the Jacobi method is known to converge very safely it suffers intrinsically from a very slow convergence for large-scale correlations (which originate mainly in the off diagonal terms of the matrix to be inverted) Press et al. (1992). This is illustrated on figure 6: there we show the maps of residuals after solving this system using a *standard* Jacobi method on simulated data with $`50`$, $`100`$, $`150`$, and $`200`$ iterations. We used the same simulation and therefore the same sky coverage. Obviously the largest structures are the slowest to converge (implying observed large scale residual patterns). As a consequence it seems very natural to help the convergence by solving the problem at larger scales. Whereas large-scale structures will not be affected by a scale change, smaller structures will converge faster.

### 4.2 Scalings

We have found that this multi-grid algorithm translates naturally in *a speed up greater than $`10`$* as compared to a standard Jacobi method. This is illustrated in figure (7) where we plotted the evolution of the 2-norm of residuals for the two methods in terms of the number of iterations in ’cycle units’. One cycle corresponds to 8 iterations at level $`k_{max}`$ for a standard Jacobi method whereas it incorporates addionally going up and down all the lowest levels in the multi-grid case. Thus the cycle timing is not exactly the same but the difference is negligible since the limiting stages are definitely the iterations performed at maximum resolution. Note the fact that the efficiency of the multi-grid method allows us to solve exactly the system up to the machine accuracy (small rebounds at the bottom of the curve) in approximately $`135\text{mn}`$ for $`(𝒩_{tod}\mathrm{8\; 10}^6,𝒩_{pix}\mathrm{8\; 10}^5)`$ on a SGI ORIGIN 2000 single processor. Since the limiting stages are the FFT’s at the higher levels, this algorithm scales as $`𝒪(𝒩_{tod}\mathrm{ln}𝒩_{tod})`$. In terms of memory storage it scales naturally as $`𝒪(𝒩_{tod})`$ since one crucial feature of this iterative method is to handle only vectors and never an explicit matrix. The scaling of this problem is formally independent of the number of pixels $`𝒩_{pix}`$.

### 4.3 The noise estimation issue and noise covariance matrix estimation

The estimation of the statistical properties of the noise in the timestream is an important issue for this kind of algorithm. Indeed, whereas till now we have assumed a perfect prior knowledge of the noise power spectrum in order to define the filters, it might not be that easy in practice since we can not separate the signal from the noise. We will therefore aim at making a joint estimation of the signal and the noise. This has been pioneered recently by Ferreira and Jaffe (2000) and implemented independently by Prunet et al. (2000). The latter implementation is rather straightforward in our particular case since it just implies reevaluating the filters after a number of iterations, given the (current) estimation of the signal map and thus of the signal data stream. Nevertheless its non-obvious convergence properties have to be studied carefully through simulations. Making use of (16) our evaluation of the noise timeline $`\widehat{n}^n`$ at the $`n^{th}`$ iteration and at level max is

$$\widehat{n}^n=dA(\widehat{y}_{k_{max}}^n+Pd).$$
(23)
Then we compute its spectrum and (re-) define the required filters. We then go through several multi-grid cycles (5 in the above demonstrated case) before re-evaluating the noise stream. Very few evaluations of the noise are needed before getting a converged power spectrum (around 2). In such an implementation, no noise priors at all are assumed. This is illustrated on one particular worked out example in the case of a 4 hours ARCHEOPS like flight (more detailed considerations will be discussed somewhere else). To reduce the number of degrees of freedom we bin the evaluated noise power spectrum using a constant logarithmic binning ($`\mathrm{\Delta }\mathrm{ln}f=0.15`$ in our case) for $`f2f_{knee}`$ and a constant linear binning ($`\mathrm{\Delta }f=0.08\mathrm{Hz}`$ in our case) for higher frequency. The figure 8 shows the genuine and evaluated noise power spectrum. The initial noise power spectrum was a realistic one $`P(f)(1+(f_{knee}/f)^\alpha )`$ to which we added some small perturbations (the two visible bumps) to test the method. Note the small bias around the telescope spin frequency at $`f_{spin}=0.05\mathrm{Hz}`$: this is illustrative of the difficulties we fundamentally face to separate signal and noise at this particular frequency through equation (16) . Naturally, this bias was not present in the case demonstrated on figure 5 where we assumed a prior knowledge of the spectrum. This possible bias forced us to work with a coarser binning ($`\mathrm{\Delta }\mathrm{ln}f=1.`$) in the $`1/f`$ part of the spectrum till the convergence is reached, i.e. we evaluate the noise power spectrum with the previously mentioned binning only at the last step. Proceeding this way, the convergence towards the correct spectrum is both fast ($`3`$ noise evaluations) and stable.

Second, the output of any map-making should contain as well an evaluation of the map noise covariance matrix $`(A^TN^1A)^1`$. Given such a fast algorithm and given an evaluation of the power spectrum, it is natural to obtain it through a Monte-Carlo algorithm fueled with various realizations of the noise timeline generated using the evaluated power spectrum. This part will be presented in a future work. However we illustrate it very briefly by a very rough $`C_{\mathrm{}}`$ determination (which is in no way an appropriate $`C_{\mathrm{}}`$ estimate). To this purpose we perform a one day TopHat like simulation including only the CMB signal plus the noise. From this data stream we obtain an “optimal” signal map as well as an evaluation of the noise power spectrum using the previously described algorithm. With the help of the anafast routine of the HEALPix package we calculate this way a rough $`C_{\mathrm{}}^{signal}`$. Using the estimated noise power spectrum we generate $`10`$ realisations of the noise and get consequently $`10`$ “optimal” noise maps. For each of them we measure as before $`C_{\mathrm{}}`$ and average them to obtain $`C_{\mathrm{}}^{noise}`$. In order to debiase the signal power spectrum recovered in this way, we substract $`C_{\mathrm{}}^{noise}`$ to $`C_{\mathrm{}}^{signal}`$. The power spectrum obtained in this manner includes as well some spatial correlations due to some light residual stripping and thus does not correspond to white noise (at least in the low $`\mathrm{}`$ part). For the aim of comparison we compute the power spectrum of the input signal measured the same way, $`C_{\mathrm{}}^{input}`$, and plotted both $`C_{\mathrm{}}^{input}`$ as well as $`C_{\mathrm{}}^{signal}C_{\mathrm{}}^{noise}`$ averaged in linear bands of constant width $`\mathrm{\Delta }\mathrm{}=80`$. The agreement is obviously very good as illustrated on figure 9. The error bars take into account both the sampling induced variance as well as the beam spreading Knox (1995). A few comments need to be made at this point. This is in no way an appropriate $`C_{\mathrm{}}`$ measurement since we are not dealing properly here with the sky coverage induced window function which triggers some spurious correlations within the $`C_{\mathrm{}}`$s. We thus do not take into account the full-structure of the map noise covariance matrix. Nevertheless, the efficiency of such a rough method is encouraging for more detailed implementation and a full handling of the noise covariance matrix.

Note finally that this is a naturally parallelised task which should therefore be feasible very quickly.

### 4.4 Application to genuine data and hypothesis

The application to genuine data could be problematic if our key hypothesis were not to be fulfilled thus we have to discuss them. Concerning the noise, we assumed that it is Gaussian in order to derive our map estimator and stationary in order to exploit the diagonality of its noise covariance matrix in Fourier space. Both hypothesis are reasonable for a nominal instrumental noise, at least partially on (sufficiently long) pieces of timeline. If not, we would have to cut the bad parts and replace them by a constrained realization of the noise in these parts. Concerning the signal, no particular assumptions are needed in the method we are presenting. At this level, we neglected as well the effect of the non perfect symmetry of the instrumental beam. This effect should be quantified for a real experiment Wu et al. (2000); Doré et al. (2000). Another technical hypothesis is the negligibility of the pixelisation noise with respect to the instrumental noise but since this is fully under control of the user it should not be a problem.

In this paper, we have presented an original fast map-making method based on a multi-grid Jacobi algorithm. It naturally entails the possibility of a joint noise/signal estimation. We propose a way to generate the noise covariance matrix and illustrate its ability on a simple $`C_{\mathrm{}}`$ estimation. The efficiency of this method has been demonstrated in the case of two coming balloon experiments, namely ARCHEOPS and TopHat but it naturally has a more general range of application. This tool should be of natural interest for Planck and this will be demonstrated somewhere else. Furthermore, due to the analogy of the formalisms, this should have some applications as well in the component separation problem.

We hope to apply it very soon on ARCHEOPS data. The FORTRAN 90 code whose results have been presented is available at http://ulysse.iap.fr/download/mapcumba/.

###### Acknowledgements.

We would like to acknowledge Eric Hivon for fruitful remarks, Lloyd Knox for stimulating discussions and information on the TopHat experiment and the developers of the HEALPix Górski et al. (1998) package which has considerably eased this work.