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

14 Pages·0.47 MB·English

by  O. Doré

A&A manuscript no. ASTRONOMY
(will be inserted by hand later) AND
Your thesaurus codes are: ASTROPHYSICS
02(12.03.1;03.13.2)
MAPCUMBA : a fast iterative multi-grid map-making
algorithm for CMB experiments
O. Dor´e1, R. Teyssier1,2, F.R. Bouchet1, D. Vibert1, and S. Prunet3
1 Institut d’Astrophysiquede Paris, 98bis, Boulevard Arago, 75014 Paris, FRANCE
2 Service d’Astrophysique,DAPNIA,Centred’E´tudes deSaclay, 91191 Gif-sur-Yvette,FRANCE
1 3 CITA, McLennan Labs, 60 St George Street, Toronto, ON M5S 3H8, CANADA
0
0 thedate of receipt and acceptance should beinserted later
2

n
a Abstract. The data analysis of current Cosmic work we are presenting focus on the first of these issues,
J
Microwave Background (CMB) experiments like namely the map-making step.

8
BOOMERanG or MAXIMA poses severe challenges Up to the most recent experiments, maps could be

1 which already stretch the limits of current (super-) made by a brute force approach amounting to solve di-
v computer capabilities, if brute force methods are used. In rectly a large linear problem by direct matrix inversion.
2 this paper we present a practical solution to the optimal Nevertheless the size of the problem, and the required
1 map making problem which can be used directly for computing power, grows as a power law of the data set
1
next generation CMB experiments like ARCHEOPS and size, and the limits of this kind of method have now been
1
0 TopHat, and can probably be extended relatively easily reached(Borrill2000).Whereasthemostefficientdevelop-
1 to the full PLANCK case. This solution is based on ment in this massive computing approach,i.e. the MAD-
0 an iterative multi-grid Jacobi algorithm which is both CAPpackage(Borrill1999)hasbeenappliedtotherecent
h/ fast and memory sparing. Indeed, if there are Ntod data BOOMERanG and MAXIMA experiments (de Bernardis
p points along the one dimensional timeline to analyse, et al. 2000;Hanany et al. 2000) some faster and less con-
- the number of operations is of O(NtodlnNtod) and the suming solutions based on iterative inversion algorithms
o
r memory requirement is O(Ntod). Timing and accuracy have now been developed and applied too (Prunet et al.
st issues have been analysed on simulated ARCHEOPS and 2000). This is in the same spirit as (Wright et al. 1996),
a TopHat data, and we discuss as well the issue of the joint and we definitely follow this latter approach.
: evaluation of the signal and noise statistical properties.
v Designgoalsare to use memory sparinglyby handling
i
X only columns or rows instead of full matrixes and to in-
crease speed by minimising the number of iterations re-
r
a Key words: methods: data analysis – cosmic microwave quired to reach the sought convergence of the iterative
background
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
1. Introduction
investigated the possibility of a joint noise (and signal)
As cosmology enters the era of “precision”, it enters si- evaluation using this algorithm.
multaneously the era of massive data sets. This has in Section 2 presents in detail the method and its imple-
turn showed the need for new data processing algorithm. mentation, while section 3 demonstrates its capabilities
Present and future CMB experiments in particular face by using simulationsoftwoon-goingballoonCMB exper-
somenewandinterestingchallenges(Bondetal.1999).If iments, ARCHEOPS1(Benoit et al. 2000) and TopHat2.
we accept the now classical point of view of a four steps The results are discussed in section 4, together with the
data analysis pipeline : i/ from time-ordered data (TOD) problem of the evaluation of the noise properties, as well
tomapsoftheskyatagivenfrequency,ii/fromfrequency as possible extensions.
maps to (among others) a temperature map, iii) from a
temperature map to its power spectrum C , iv/ from the
ℓ
power spectrum to cosmological parameters and charac-
teristics of the primordialspectrum of fluctuation, the ul- 2. The method
timate quantities to be measured in a given model. The
Send offprint requests to: O. Dor´e 1 http://www-crtbt.polycnrs-gre.fr/archeops/
Correspondence to: [email protected] 2 http://topweb.gsfc.nasa.gov/
2 Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code
2.1. Optimal map making the least square solution which was used to analyse the
“COBE” data (Jansen and Gulkis 1992),
The relation between the sky map we seek and the ob-
serveddatastreammaybecastasalinearalgebrasystem W =(cid:2)ATN−1A(cid:3)−1ATN−1. (8)
(Wrightetal.1996;Tegmark1997).Let and indicesde-
t p
note quantities in the temporal and spatial domains, and In this paper we will actually deal only with this estima-
group as a data vector, d , and a noise vector the tempo- tor. Nevertheless as a next iteration in the analysis pro-
t
ralstreamofcollecteddataandthedetectornoisestream, cess, we could incorporate various theoretical priors by
both of dimension N . We then have explicitingP(x).Forexample,itisoftenassumedaGaus-
tod
sian prior for the theory, i.e. P(x)∝exp(−xTpCp−p1′xp′/2)
dt =Atpxp+nt, (1) where Cpp′ = hxpxTp′i is the signal covariance matrix.
In that case the particular solution turns out to be the
where A x is the signal vector given by the observation
tp p
Wienerfilteringsolution(Zaroubietal.1995;Bouchetand
ofthe unknownpixelisedskymap,x ,whichhasbeenar-
p
Gispert1996;TegmarkandEfstathiou1996;Bouchetand
ranged as a vector of dimension N . The N ×N
pix tod pix
Gispert 1998):
“observation” matrix A therefore encompasses the scan-
ning strategy and the beam pattern of the detector. W =(cid:2)C−1+ATN−1A(cid:3)−1ATN−1 . (9)
Inthefollowing,werestricttothecasewhenthebeam
pattern is symmetrical. We can therefore take xp to be a But this solution may always be obtained by a further
map of the sky which has already been convolved with (Wiener) filtering of the COBE solution, and we do not
the beam pattern, and A only describes how this sky is consider it further.
beingscanned.Forthetotalpowermeasurement(i.e.non- The prior-less solution demonstrates that as long as
differential measurement) we are interested in here, the the (Gaussian) instrumental noise is not white, a sim-
observation matrix A then has a single non-zero element ple averaging (co-addition) of all the data points corre-
per row, which can be set to one if d and x are expressed sponding to a given sky pixel is not optimal. If the noise
in the same units. The model of the measurement is then exhibits some temporal correlations, as induced for in-
quite transparent:eachtemporaldatum is the sum of the stance by a low-frequency 1/f behaviorof the noise spec-
pixel value to which the detector is pointing plus the de- trum which prevails in most CMB experiments, one has
tector noise. to take into account the full time correlation structure of
Themap-makingstepthenamountstobestsolveforx the noise. Additionally this expression demonstrates that
given d (and some properties of the noise). We shall seek even if the noise has a simple time structure, the scan-
a linear estimator of x , ning strategy generically induces a non-trivial correlation
p
matrix (cid:2)ATN−1A(cid:3)−1 of the noise map.
xˆ =W d . (2)
p pt t Eveniftheproblemiswellposedformally,aquicklook
attheordersofmagnitudeshowsthattheactualfindingof
To motivate a particularchoiceof theN ×N matrix
pix tod
asolutionisnontrivialtask.Indeedabruteforcemethod
Wing,athBeaoypestiiamnaalpsporluotaicohnistocotnhviseninievnetr.sIinondepedrowbleeamrewsheeickh- aiming at inverting the full matrix (cid:2)ATN−1A(cid:3)−1, an op-
eration scaling as O(N3 ), is already hardly tractable
maximises the probability of a deduced set of theory pa- pix
for present long duration balloon flights as MAXIMA,
rameters (here the map x ) given our data (d ) by max-
p t
BOOMERanG, ARCHEOPS or TopHat where N ∼
imising P(x|d). Bayes’ theorem simply states that tod
106 and N ∼ 105. It appears totally impractical for
pix
P(d|x)P(x) PLANCKsinceforasingledetector(amid10s)N ∼109
P(x|d)= . (3) tod
P(d) and Npix ∼107!
One possibility may be to take advantage of specific
If we do not assume any theoretical prior, then x follows
scanning strategies, and actually solve the inverse of the
a uniform distribution as well as d. Therefore,
convolution problem as detailed in (Wandelt and Gorski
2000). This amounts to deduce the map coefficients a
P(x|d) ∝ P(d|x). (4) lm
in the spherical harmonic basis through a rotation of a
If we further assume that the noise is Gaussian, we can Fourierdecompositionoftheobserveddata.Themapwill
write thenbe asimplevisualisationdevice,whilethealm would
be ready to use directly for a component separation (as
P(x|d) ∝ exp(−nTt Nt−t′1nt′/2) (5) in (Bouchet and Gispert 1996; Tegmark and Efstathiou
∝ exp(−(d−Ax)Tt Nt−t′1(d−Ax)t′/2) (6) 1996; Bouchet and Gispert 1998)) and the CMB power
∝ exp(−χ2/2) (7) spectrumestimate.Whilepotentiallyveryinteresting,this
approach will not be generally applicable (at least effi-
whereNtt′ =<nnT >tt′ isthenoisecovariancematrix.In ciently),andwenowturntoapractical(general)solution
this particular case, maximising P(x|y) amounts to find of equation (8) by iterative means.
Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 3
2.2. Practical implementation for large data sets which will be non-zero in practice, as soon as one works
withfinitespatialresolution.Wecallthisresidualthepix-
We solve the map-making problem by adapting to our
elisation noise. Since we assume here that the instrumen-
particular case the general“multi-grid method” (Press et
tal beam is symmetric, the sky map is considered as the
al.1992).Multi-gridmethodsarecommonlyusedtospeed
true sky convolved by, say, a Gaussian beam of angular
up the convergenceof a traditionalrelaxationmethod (in
diameter ∆θ . This introduces a low-pass spatialfilter in
B
our case the Jacobi method, as in (Prunet et al. 2000))
theproblem.Inotherwords,astheresolutionincrease,the
definedatresolutionℓ (seebelow).Asetofrecursively
max pixelisation noise should decrease towards zero. We have
defined coarser grids (ℓ < ℓ ) are used as temporary
max estimated that the order of magnitude of the pixelisation
computationalspace,inordertoincreasethe convergence
noise can be approximated by
rate of the relaxation process. To be fully profitable, this
algorithm implies for each resolution both a rebinning in ∆θk
space (resolution change) and in time (resampling). kpkk≃kx∞k∆θ (14)
B
In this paper, we use the HEALPix pixelisation of the
The norms used in the above formula can be either the
sphere (Go´rski et al. 1998). In this scheme, the sphere
maximum over the time line (a very strong constraint) or
is covered by 12 basic quadrilaterals, further divided re-
thevarianceoverthetimeline(aweakerconstraint).Since
cursively into pixels of equal area. The map resolution
the pixelisation noise is strongly correlated with the sky
is labeled by N : the number of pixels along the side
side
signal,pointsourcesorGalaxycrossingsarepotentialcan-
of one basic quadrilateral. Hence, N = 1 means that
side
didates forlargeandlocalisedburstsofpixelisationnoise.
the sphere is coveredby 12 large pixels only. The number
The correct working resolution k is set by requiring
of pixels is given by N = 12N2 . N = 256 corre- max
pix side side that the pixelisation noise remains low compared to the
spondstoapixelsizeof13.7arcmin.Forpracticalreasons,
actual instrumental noise, or equivalently
we need to define the resolution k of a HEALPix map as
∆θ
B
Nside =2k (10) ∆θkmax ≤ S/N (15)
The“nested”pixelnumberingschemeofHEALPix(Go´rski Most of the CMB experiments are noise dominated along
etal.1998)allowsaneasyimplementationofthe coarsen- the time line, constraining the effective map resolutionto
ing (k → k−1) and refining (k → k+1) operators that be of the order of the instrumental beam or even larger.
we use intensively in our multigrid scheme. Note however that the pixelisation noise is strongly non
Let us now get into the details of our implementation Gaussian (point sources or Galaxy crossings) and can be
and discuss successively the exact system we solve, the alwaysconsideredas a potential source of residualstripes
way we solve it and the actual steps of the multi-grid al- in the final maps.
gorithm.
2.2.2. The basic relaxation method: an approximate
2.2.1. Determining the working resolution kmax Jacobi solver
Weaimatsolvingfortheoptimalmapxˆkatagivenspatial Insteadofsolving forxˆ we performthe changeofvariable
resolution k using
yˆ=xˆ−P d (16)
ATN−1A xˆ =ATN−1 d, (11)
k k k k
and solve instead for yˆwhich obeys
where A is the “observation” operator (from spatial to
temporalkdomain) and AT is the “projection” operator PN−1A yˆ=PN−1(d−APd) or Myˆ=b (17)
k
(from temporal to spatial domain). In a noise-free exper-
where we have multiplied each side of equation (11) by
iment, the optimal map would be straightforwardlygiven (ATA)−1, the pixel hit counter. From now on, we also
by the co-added map (introducing the “co-addition” op-
assume that the noise in the timestream is stationary
erator P )
k and that its covariance matrix is normalised so that
xˆ =P d≡(ATA )−1AT d (12) diag Nt−t′1 = I. The previous change of variable allows
k k k k k us to subtractthe sky signalfromthe data:since we have
The time line is givenby d=A∞x∞+n where x∞ is the chosenaresolutionhighenoughtoneglectthepixelisation
sky map at “infinite” resolution (k = +∞ in our nota- noise, we have indeed APd ≃ A∞x∞+APn and, conse-
tions).In orderto check the accuracyofthis trivialnoise- quently,d−APd≃n−APn.Themapmakingconsistsin
free map making, it is natural to compute the residual in two step: first compute a simple co-added map from the
the time domain with n=0 time line, and second, solve equation (17) for the stripes
mapyˆ.Thefinaloptimalmapisobtainedbyaddingthese
pk =Akxˆk−d=Akxˆk−A∞x∞ (13) two maps.
4 Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code
It is worth mentioning that the stripes map is com- Since both A and P are norm-preserving operators, the
pletely independent of the sky signal, as soon as the pix- norm of the increment ∆yˆn = yˆn − yˆn−1 between step
elisationnoisecanbeignored.Itdependsonlyonthescan- n and n + 1 decreases as k∆yˆn+1k ≤ W(f )k∆yˆnk,
min
ningstrategythroughthematrixAandonthenoisechar- where f is a minimal frequency in the problem. Since
min
acteristicsthroughthe matrix N.Evenif inprinciple this W(f ) < 1, we see that the approximate Jacobi relax-
min
changeofvariableisirrelevantsinceitdoesnotchangethe ation scheme will converge in every case, which is good
matrix to be inverted, it does in practice since d−APd news. On the other hand, since W(f ) ≃ 1, the ac-
min
is free from the fast variationsof d (up to the pixelisation tual convergence rate of the scheme is likely to be very,
noise),ase.g. theGalaxycrossings,whicharenumerically very slow, which is bad news (cf. figure 6 for a graphi-
damaging to the quality of the inversion. cal illustration) . The fact that this algorithm is robust,
Tosolveforequation17,wefollowtheideasof(Prunet but dramatically slow is a well-knownproperty of the Ja-
et al. 2000) and apply an approximate Jacobi relaxation cobi method. The multi-grid method is also well known
scheme. The Jacobi method consists in the following re- to solve this convergence speed problem. Note that if the
laxation step to solve for our problem convergenceis reached,the solution we get is the optimal
solution,i.e. similartotheonethatwouldbe obtainedby
yˆn+1 =Ryˆn+D−1b and yˆ0 =0 (18)
a full matrix inversion.
where D is the diagonal part of the matrix M and R is
the residualmatrix R=I−D−1M.Computing the diag- 2.2.3. Multigrid relaxation
onal elements of M = PN−1A is rather prohibitive. The
idea of Prunet et al. (2000) is to approximate D ≃ I by The multi-grid method for solving linear problems is de-
neglecting the off-diagonalelements of N−1. The residual scribed in great details in (Press et al. 1992). At each
matrix then simplifies greatly and writes relaxation step at level k = kmax, our target resolution,
we can define the error en = yˆn − yˆ and the residual
R=P(I −N−1)A (19) rn =M yˆn−b . Both arekrelatedk throkugh
k k k k
The approximate Jacobi relaxation step is therefore de-
M en =rn (22)
fined as k k k
yˆn+1 =Ryˆn+b and yˆ0 =0 (20) Ifweareabletosolveexactlyforequation(22),theoverall
problem is solved. The idea of the multi-grid algorithmis
One clearly sees that if this iterative scheme converges, to solve approximately for equation (22) using a coarser
it is towards the full solution of equation 17. To perform “grid” at resolution k −1, where the relaxation scheme
thesesuccessivesteps,itisextremelyfruitfultoremember should converge faster. We thus need to define a fine-to-
theassumedstationarityofthenoise.Indeedwhereasthis coarse operator, in order to define the new problem on
assumption implies a circulant noise covariance matrix in the coarse grid and solve for it. We also need a coarse-to-
realspace,ittranslatesinFourierspaceinthediagonality fine operator in order to inject the solution onto the fine
of the noise covariance matrix. This is naturally another grid. The approximation to the error en is finally added
k
formulationoftheconvolutiontheorem,sinceastationary to the solution. The coarse grid solver is usually applied
matrix acts on a vector as a convolution, and a diagonal afterafewiterationsonthefinelevelhavebeenperformed
matrixactsasasimplevectorproduct,thusaconvolution (inpracticeweperform3to5iterations).Naturally,since
in real space is translated as a product in Fourier space. the solution to the problem on the coarse level relies also
The point is that the manipulation of the matrix N−1 is on the same relaxation scheme, it can be itself acceler-
considerably lighter and will be henceforth performed in ated using an even coarser grid. This naturally leads to a
Fourier space. Applying the matrix R to a map reduces recursive approach of the problem.
then to the following operations in order We defined our fine-to-coarse operator to be an aver-
aging of the values of the 4 fine pixels contained in each
1. “observe” the previous stripes map yˆn
coarse pixel. The coarse-to-fine operator is just a straight
2. Fourier transform the resulting data stream
injection of the value of the coarse pixel into the 4 fine
3. apply the low-pass filter W =I−N−1
pixels. The most important issue is the temporal rebin-
4. inverse Fourier transform back to real space
ning of the data stream since the speed of the iterative
5. co-add the resulting data stream into the map yˆn+1.
scheme at a given level is actually set by the two Fourier
transforms.Weperformedthatresamplingbysimplytak-
Assuming that the normalised noise power spectrum
ing eachtime we go up one level one point out of two. At
α
canbe approximatedbyP(f)=1+(f0/f) ,the low-pass thelowerresolutions,thisreductionissuchthattheitera-
filter associated to each relaxation step is given by
tion cost is negligible when comparedto that at higher k;
fα it allows fast enough iterations that full convergence may
W(f)= 0 (21)
fα+fα be reached. In practice we choose a minimal level k = 3
0
Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 5
defined by N = 8 and iterate a few hundred times to Weintroduced5distinctlevelsofresolutiondefinedby
side
reach exact convergence. their N parameter in the HEALPix package (Go´rski
side
Finallythenavigationthroughlevelsallowsseveralop- et al. 1998). The higher resolution level is imposed by
tions to be taken. Either we go up and down through all the pixelisation noise level requirement (section 2.2.3) to
′
the levels successively (the so-called “V-cycle”) or we fol- N = 256 (pixel size ∼ 13.7) whereas the lower one is
side
low more intricate paths (e.g. the “W-cycle” where at a N = 8 (pixel size ∼ 7.3o). Therefore these two config-
side
given level we go down and up all the lower levels before urations each offer an interesting test since they differ by
goingupandvice-versa).Sinceitturnsoutthatlevelsare the sky coverage and the noise power spectrum. We iter-
relatively disconnected in the sense that the scales well ate 3 times at each level except at the lowest one where
solved at a given resolution are only slightly affected by we iterate 100 times.
thesolutionatahigherlevel,the“V-cycle”isappropriate
and is the solutionwe adopt in the two following configu-
3.2. Results
rations.
The algorithm is as efficient in both situations. Whereas
forARCHEOPSwhosetimelineislongerduetothehigher
3. Practical application to ARCHEOPS and
sampling frequency it took 2.25 hours on a SGI ORIGIN
TopHat experiments
2000 single processor, it took around 1.5 hours for the
We now aim at demonstrating the capabilities of this al- TopHat daily data stream.
gorithmwith simulated data whose characteristicsare in- In figures 1 and 2 we depict from top to bottom and
spired by those of the ARCHEOPS and TopHat experi- from left to right the initial co-added data, the recon-
ments. structed noise map, the hit map, i.e. the number of hit
per pixels atthe highestworkedoutresolution,the initial
3.1. Simulation signalmapaswellasthereconstructedsignalmapandthe
error map. We see that the destriping is excellent in both
The ARCHEOPS and TopHat experiments are very sim- situations and the signal maps recovered only contain an
ilar with respect to their scanning strategy. Indeed both
apparently isotropic noise. We note the natural correla-
use a telescope that simply spins at a constant rate (re-
tion between the error map and the hit map. Finally, we
spectively3and4RPM)aboutitsverticalaxis.Thusdue stress that no previous filtering at all has been applied.
to Earth rotation the sky coverage of both is performed
Figure 3 shows how the noise map is reconstructed at
through scan circles whose axis slowly drifts on the sky.
variousscales.Thisisamereillustrationofourmulti-grid
Nevertheless, because of the different launch points (re-
work.
spectively on the Arctic circle (Kiruna, Sweden) and in
Antarctica (McMurdo)) and their different telescope axis
elevation(respectively∼41o or12o)theydo nothavethe 3.3. Tests
same sky coverage.
At this point, some tests are required to cautiously vali-
Otherwise the two experiments do not use the same
date our method. First, as was stated below, as soon as
bolometers technology, neither the same bands nor have
the iterative algorithm has converged, the solution is by
the same number of bolometers. But even if we try to be
construction the optimal solution, similar to the one that
fairlyrealistic,ourgoalthoughisnottocomparetheir re-
would be obtained by the full matrix inversion. As a cri-
spective performances but rather to demonstrate two ap-
teriumforconvergencewerequiredthe2-normofresiduals
plications of our technique in different settings. We then
to be of the order of the machine accuracy.
simulate for each a data stream of ∼ 24 hrs duration
WeinitiallyhaveaGaussianrandomnoisestreamfully
with respectively a sampling frequency of 171 and 64 Hz.
characterisedby its powerspectrum.Thereforean impor-
The TODs contain realistic CMB and Galactic signal for
tant test is to check wether the deduced noise stream (by
a band at 143 GHz. Note that this is a one day data
“observing” the stripe map, see § 4.3 for further details
streamfor TopHat (out of 10 expected) and that this fre-
definition)is Gaussianandhasthe samepowerspectrum.
quency choice is more appropriate for ARCHEOPS than
Onfigure4weensurethattheevaluatednoisetimestream
for TopHat (whose equivalent band is around 156 GHz),
is indeed gaussian.As depicted onfigure 5,where we plot
but this is realistic enoughfor our purpose.We generated
in the ARCHEOPScaseboth the analyticalexpressionof
a Gaussian noise time stream with the following power
the spectrumaccordingtowhichwegeneratethetimeline
spectrum P(f) ∝ (1+(f /f)α)−1 with f = 0.24
knee knee and its recoveredform, the agreementis excellent. We re-
and 1 Hz, and α = 1.68 and 1. The noise amplitude per
callthatweassumedatthebeginningaperfectknowledge
channelischoosensothatitcorrespondsforARCHEOPS
of the noise inorder to define the filters.This is naturally
(24 hours)andTopHat3 (10days offlight)to 30/8µK on
unrealisticbuttheissueofnoiseevaluationisdiscussedin
′ ′ ′
average per 20 pixel, with a beam FWHM of 10 / 20.
section 4.3 below. We plotted as well the probability dis-
3 Lloyd Knox privatecommunication tribution function (PDF) of the final error map, i.e. the
6 Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code
Fig.1. Simulated ARCHEOPS Kiruna flight : From top to bottom and from left to right, the co-added map and the
inputGalactic+ CMBsignalmap,the reconstructednoise(stripes)andsignalmaps,the hitcountandtheerrormap
(Arbitrary unit). The fact thatthe coverageis notfully symmetricalis due to the fact that we consideredslightly less
′
than 24hr. Mollweide projection with pixels of 13.7 (HEALPix N =256). Arbitrary units.
side
Fig.2. Simulated TopHat one day flight : From top to bottom and from left to right, the co-added map and the
input CMB + Galaxy signal map, the reconstructed noise (or stripes) and signal map, the hit count and error map.
The factthat the coverageis notfully symmetricalis due to the fact that we consider only 18.2hrof flight.Gnomonic
′
projection with pixel of 13.7 (HEALPix N = 256). Note that the slight visible stripping is correlated to the
side
incomplete rotation pattern. Arbitrary units.
Fig.3. Multi-grid noise map recovery : In this plot we show in the ARCHEOPS Kiruna case how the noise map is
reconstructed at various levels, corresponding respectively to N =8,16,32,64.
side
recovered noisy signal map minus the input signal map 4.2. Scalings
(figure 4). This PDF is well fitted by a gaussian whose
parameters are given in figure 4. The PDF of the error We have found that this multi-grid algorithm translates
map displays some non-gaussianwings. Let us recall here naturally in a speed up greater than 10 as compared to a
that this is no surprise here because of the non-uniform standard Jacobi method. This is illustrated in figure (7)
sky coverage as well as the slight residual stripping, both where we plotted the evolution of the 2-norm of residu-
due to the non-ideal scanning strategy, i.e. that produces als for the two methods in terms of the number of iter-
a non-uniform white noise in pixel space. ations in ’cycle units’. One cycle corresponds to 8 itera-
tions at levelk for a standardJacobimethod whereas
Another particularly important test consists in check- max
it incorporatesaddionally going up and down all the low-
ingtheabsenceofcorrelationbetweenthe recoverednoise
est levels in the multi-grid case. Thus the cycle timing
map and the initial signal map. We could not find any
is not exactly the same but the difference is negligible
which is no surprise since we are iterating on a noise map
since the limiting stages are definitely the iterations per-
(see § 2.2.2) which does not contain any signal up to the
formedatmaximumresolution.Notethefactthattheeffi-
pixelisationnoise,thatisensuredtobenegligiblewithre-
ciency of the multi-gridmethodallows us to solve exactly
gardstotheinstrumentalnoisebychoiceoftheresolution
the system up to the machine accuracy (small rebounds
k (see § 2.2.1).
max
at the bottom of the curve) in approximately 135mn for
(N ∼8 106,N ∼8 105) on a SGI ORIGIN 2000 sin-
tod pix
gle processor. Since the limiting stages are the FFT’s at
4. Discussion
the higher levels,this algorithmscalesas O(N lnN ).
tod tod
4.1. Why is such an algorithm so efficient ? IntermsofmemorystorageitscalesnaturallyasO(N )
tod
sinceonecrucialfeatureofthisiterativemethodistohan-
The efficiency of such a method can be intuitively un-
dle only vectorsandnever anexplicitmatrix.The scaling
derstood. Indeed, although the Jacobi method is known
of this problem is formally independent of the number of
to converge very safely it suffers intrinsically from a very
pixels N .
pix
slow convergence for large-scale correlations (which origi-
nate mainly in the off diagonal terms of the matrix to be
inverted) (Press et al. 1992). This is illustrated on figure 4.3. The noise estimation issue and noise covariance
6: there we show the maps of residuals after solving this matrix estimation
systemusingastandard Jacobimethodonsimulateddata
with 50, 100, 150, and 200 iterations. We used the same The estimationofthe statisticalpropertiesofthe noise in
simulationandthereforethesameskycoverage.Obviously the timestreamisanimportantissueforthiskindofalgo-
the largest structures are the slowest to converge (imply- rithm.Indeed,whereastillnowwehaveassumedaperfect
ing observed large scale residual patterns). As a conse- prior knowledge of the noise power spectrum in order to
quence it seems very natural to help the convergence by definethefilters,itmightnotbethateasyinpracticesince
solving the problem at larger scales. Whereas large-scale we can not separate the signal from the noise. We will
structures will not be affected by a scale change, smaller therefore aim at making a joint estimation of the signal
structures will converge faster. and the noise. This has been pioneered recently by (Fer-
Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 7
Fig.6. Residual map after 50, 100, 150 and 200 iterations of a standard Jacobi method. This has been performed on
simulated data for one bolometer with a nominal noise level.The sky coverageis that of ARCHEOPScoming Kiruna
flight. The residual large scale patterns illustrate the difficulties the standard Jacobi method faces to converge. The
stripes free area are the ones of scan crossing (see the hit map in figure 1).
Fig.4. In the TopHat case we plot from top to bottom Fig.5. Recovery of the noise power spectrum in the
therecoveredprobabilitydistributionfunctionofthenoise Archeopscase(top:linearx-axis,bottom:logx-axis):The
stream evaluated along the timeline as well as the error redthindashedlineshowsthe initialanalyticnoisepower
map PDF. In this two cases a fit to a gaussian has been spectrumusedtogeneratethe inputnoisestreamandthe
performedwhoseparametersarewritteninsidethefigures. bluethicklinedenotestherecoveredoneafter6iterations.
No significant departure from gaussianity are detected. Therecoveredonehasbeenbinnedasdescribedinsection
The arbitrary units are the same as the ones used for the 4.3andbothhavebeennormalisedso thatthe white high
previously shown maps. frequency noise spectrum is equal to one. The agreement
is obviously very good. No apparent bias is visible. Note
that a perfect noise power spectrum prior knowledge has
been used in this application.
reira and Jaffe 2000) and implemented independently by
(Prunet et al. 2000). The latter implementation is rather
straightforwardin our particularcase since it just implies
reevaluating the filters after a number of iterations, given simulations.Makinguseof(16)ourevaluationofthenoise
the (current) estimation of the signal map and thus of timeline nˆn at the nth iteration and at level max is
the signal data stream. Nevertheless its non-obvious con-
vergence properties have to be studied carefully through nˆn =d−A(yˆn +Pd). (23)
kmax
8 Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code
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 conver-
gence towards the correct spectrum is both fast (3 noise
evaluations) and stable.
Second,the outputofanymap-makingshouldcontain
as well an evaluation of the map noise covariance matrix
(ATN−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
realizationsofthenoisetimelinegeneratedusingthe eval-
uated power spectrum. This part will be presented in a
future work. However we illustrate it very briefly by a
very rough C determination (which is in no way an ap-
ℓ
propriateC estimate). To this purpose we perform a one
ℓ
day TopHat like simulation including only the CMB sig-
Fig.7. The evolution of the 2-norm of residuals with
nal plus the noise. From this data stream we obtain an
the number of iterations at level max. Whereas the blue
“optimal”signalmapaswellasanevaluationofthe noise
dashed line is standard iterative Jacobi, the solid red line
powerspectrumusingthepreviouslydescribedalgorithm.
is the iterative multi-grid Jacobi method. A full multi-
With the help of the anafast routine of the HEALPix
grid cycle incorporates 3 iterations at level max before package we calculate this way a rough Csignal. Using the
ℓ
going down and up all the levels. The sharp jumps corre-
estimated noise power spectrum we generate 10 realisa-
sponds to the moment when the scheme reach againlevel
tionsofthenoiseandgetconsequently10“optimal”noise
max and thus take the full benefits fromthe resolutionat
maps.ForeachofthemwemeasureasbeforeC andaver-
ℓ
lowerlevels.Notethattheverysharpfeaturesafter∼100 age them to obtain Cnoise. In order to debiase the signal
ℓ
iterationsareduetothefactwereachedthemachineaccu- powerspectrumrecoveredinthisway,wesubstractCnoise
ℓ
racy which is almost out of reach for a standard iterative to Csignal. The power spectrum obtained in this manner
method. ℓ
includesaswellsomespatialcorrelationsduetosomelight
residual stripping and thus does not correspond to white
noise (at least in the low ℓ part). For the aim of compar-
Then we compute its spectrum and (re-) define the re- ison we compute the power spectrum of the input signal
quired filters. We then go through several multi-grid measured the same way, Cinput, and plotted both Cinput
ℓ ℓ
cycles (5 in the above demonstrated case) before re- as well as Csignal − Cnoise averaged in linear bands of
evaluating the noise stream. Very few evaluations of the ℓ ℓ
constantwidth ∆ℓ=80.The agreementis obviouslyvery
noise are needed before getting a converged power spec-
good as illustrated on figure 9. The error bars take into
trum (around 2). In such an implementation, no noise
accountboththesamplinginducedvarianceaswellasthe
priors at all are assumed. This is illustrated on one par-
beamspreading(Knox1995).Afewcommentsneedtobe
ticular worked out example in the case of a 4 hours
made at this point. This is in no way an appropriate C
ℓ
ARCHEOPS like flight (more detailed considerations will
measurement since we are not dealing properly here with
be discussed somewhere else). To reduce the number of
the sky coverage induced window function which triggers
degreesoffreedomwebintheevaluatednoisepowerspec-
some spurious correlations within the C s. We thus do
ℓ
trum using a constant logarithmic binning (∆lnf =0.15
not take into account the full-structure of the map noise
in our case) for f ≤ 2 f and a constant linear bin-
knee covariance matrix. Nevertheless, the efficiency of such a
ning (∆ f = 0.08 Hz in our case) for higher frequency.
roughmethodisencouragingformoredetailedimplemen-
Thefigure8showsthe genuineandevaluatednoisepower
tation and a full handling of the noise covariance matrix.
spectrum. The initial noise power spectrum was a real-
Note finally that this is a naturally parallelised task
istic one P(f) ∝ (1+(f /f)α) to which we added
knee
which should therefore be feasible very quickly.
some small perturbations (the two visible bumps) to test
themethod.Notethesmallbiasaroundthetelescopespin
frequencyatfspin =0.05Hz:thisisillustrativeofthediffi- 4.4. Application to genuine data and hypothesis
culties wefundamentallyface toseparatesignalandnoise
atthisparticularfrequencythroughequation(16).Natu- The application to genuine data could be problematic if
rally,thisbiaswasnotpresentinthecasedemonstratedon our key hypothesis were not to be fulfilled thus we have
figure 5 where we assumed a prior knowledge of the spec- to discuss them. Concerning the noise, we assumed that
trum. This possible bias forced us to work with a coarser it is Gaussian in order to derive our map estimator and
binning (∆lnf = 1.) in the 1/f part of the spectrum stationary in order to exploit the diagonality of its noise
Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 9
Fig.9. In the case of a one day flight of the TopHat ex-
periment we perform a very approximate evaluation of
the recoveredsignal band powers (black) which has to be
compared to the input signal power spectrum (red line).
Bothhavebeenmeasuredthesamewayusingtheanafast
routine of the HEALPix package. These band powers
(∆ℓ=80) have been performed using a fast Monte-Carlo
evaluationofCnoise (dashedblueline)whichdoesnotcor-
ℓ
respond exactly to white noise since there remains some
spatial correlations.This constitutes in no way an appro-
priateC measurementbutisanencouragingsteptowards
ℓ
a full Monte-Carlo treatment.
In this paper,we have presented an originalfast map-
making method based on a multi-grid Jacobi algorithm.
Fig.8. Evaluation of the noise power spectrum (top: lin- It naturally entails the possibility of a joint noise/signal
ear x-axis, bottom: log x-axis): the red thin dashed line estimation. We propose a way to generate the noise co-
shows the initial noise power spectrum obtained from the variance matrix and illustrate its ability on a simple C
ℓ
input noise streamand the blue thick line denotes the re- estimation.Theefficiencyofthismethodhasbeendemon-
covered one after 5 iterations. Both have been smoothed strated in the case of two coming balloon experiments,
and normalised so that the white high frequency noise namely ARCHEOPS and TopHat but it naturally has a
spectrumisequaltoone.Theagreementisobviouslyvery more general range of application. This tool should be of
good. Note that no noise priors at all have been used in natural interest for Planck and this will be demonstrated
this evaluation. somewhere else. Furthermore, due to the analogy of the
formalisms, this should have some applications as well in
the component separation problem.
covariance matrix in Fourier space. Both hypothesis are
We hope to apply it very soon on
reasonable for a nominal instrumental noise, at least par-
ARCHEOPS data. The FORTRAN 90 code whose
tially on (sufficiently long) pieces of timeline. If not, we
results have been presented is available at
would have to cut the bad parts and replace them by a
http://ulysse.iap.fr/download/mapcumba/.
constrained realization of the noise in these parts. Con-
cerning the signal, no particular assumptions are needed
Acknowledgements. Wewould like to acknowledge Eric Hivon
in the method we are presenting. At this level, we ne-
for fruitful remarks, Lloyd Knox for stimulating discussions
glected as well the effect of the non perfect symmetry of and information on the TopHat experiment and the develop-
the instrumental beam. This effect should be quantified ers of the HEALPix (G´orski et al. 1998) package which has
for a real experiment (Wu et al. 2000; Dor´e et al. 2000). considerably eased thiswork.
Anothertechnicalhypothesisisthenegligibilityofthepix-
elisation noise with respect to the instrumental noise but
References
since this is fully under control of the user it should not
be a problem. Benoit A.et al., submitted to A&A2000
10 Dor´e, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code
Bond J.R., Crittenden R.G., R.G., Andrew H., Jaffe A.H.,
Knox L., Computing in Science and Engineering (1999),
astro-ph/9903166
Borrill J.,‘MADCAP: The Microwave Anisotropy Dataset
Computational Analysis Package’ Julian Borrill
in ”Proceeedings of the 5th European SGI/Cray
MPP Workshop” (1999), astro-ph/9911389, see also
http://cfpa.berkeley.edu∼borrill/cmb/madcap.html
Borrill J., oral contribution totheMPA/ESO/MPE Joint As-
tronomy Conference “ Mining thesky”, Garching 2000
Bunn E.F., Hoffman Y., Silk J., ApJ, 464, 1 (1996)
BouchetF.R.,Gispert,AIPConferenceProceedings348,‘The
mm/sub-mm Foregrounds and Future CMB Space Mis-
sions”, Ed. E. Dwek (1996)
Bouchet F. R., Gispert, R., 1999, New Astron., 4, 443B
deBernardis P. et al., Nature (London) 404, 995 (2000)
Dor´eO.,BouchetF.R.,TeyssierR.,VibertD.inproceedingsto
the MPA/ESO/MPE Joint Astronomy Conference “ Min-
ing thesky”, Garching 2000
Ferreira P.G. and Jaffe A.H.,MNRAS,312, 89 (2000)
G´orski E.K., Hivon E., Wandelt B.D. in proceed-
ings of the MPA/ESO Garching Conference 1998,
eds Banday A.J., Sheth K. and L. Da Costa and
http://www.eso.org∼kgorski/healpix/
Hanany S. et al., Astrophys. J. Lett. submitted (2000), astro-
ph/0005123
JansenD.J.andGulkisS.,“MappingtheSkywiththeCOBE-
DMR”, in “The Infrared and Submillimeter Sky after
COBE”,edsM.SignoreandC.Dupraz(Dodrecht:Kluwer)
Knox L., Phys. Rev.D,52, 4307, 1995
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.,
1992, Numerical Recipes, 2nd ed., Cambridge University
Press, Cambridge.
Prunet, S. and Netterfield, C. B. and Hivon, E. and Crill, B.
P.,”Iterativemap-makingforscanningexperiments”,tobe
publishedintheProceedingsoftheXXXVthRencontresde
Moriond, Energy Densities in theUniverse,Editions Fron-
tieres (2000), astro-ph/0006052
Tegmark M., ApJLett., 480, L87-L90 (1997)
Tegmark, M. and Efstathiou, G., MNRAS,281, 1297 (1996)
Wandelt,B. D.and Gorski, K.M., preprint astro-ph/0008227
Wu J.H.P. et al., APJs in press, astro-ph/0007212
Wright E.L., Hinshaw G., Bennett C.L., ApJ Lett., 458,53
(1996)
Wright E.L., proceeding of the IAS CMB Workshop, astro-
ph/9612006
ZaroubiS.,HoffmanY.,FisherK.B.,LahavO.,ApJ,449,446
(1995)