We present a multi-level solver for drawing constrained Gaussian realizations or finding the maximum likelihood estimate of the CMB sky, given noisy sky maps with partial sky coverage. The method converges substantially faster than existing Conjugate
Gradient (CG) methods for the same problem. For instance, for the 143 GHz Planck frequency channel, only 3 multi-level W-cycles result in an absolute error smaller than 1 microKelvin in any pixel. Using 16 CPU cores, this translates to a computational expense of 6 minutes wall time per realization, plus 8 minutes wall time for a power spectrum-dependent precomputation. Each additional W-cycle reduces the error by more than an order of magnitude, at an additional computational cost of 2 minutes. For comparison, we have never been able to achieve similar absolute convergence with conventional CG methods for this high signal-to-noise data set, even after thousands of CG iterations and employing expensive preconditioners. The solver is part of the Commander 2 code, which is available with an open source license at http://commander.bitbucket.org/.
We describe and implement an exact, flexible, and computationally efficient algorithm for joint component separation and CMB power spectrum estimation, building on a Gibbs sampling framework. Two essential new features are 1) conditional sampling of
foreground spectral parameters, and 2) joint sampling of all amplitude-type degrees of freedom (e.g., CMB, foreground pixel amplitudes, and global template amplitudes) given spectral parameters. Given a parametric model of the foreground signals, we estimate efficiently and accurately the exact joint foreground-CMB posterior distribution, and therefore all marginal distributions such as the CMB power spectrum or foreground spectral index posteriors. The main limitation of the current implementation is the requirement of identical beam responses at all frequencies, which restricts the analysis to the lowest resolution of a given experiment. We outline a future generalization to multi-resolution observations. To verify the method, we analyse simple models and compare the results to analytical predictions. We then analyze a realistic simulation with properties similar to the 3-yr WMAP data, downgraded to a common resolution of 3 degree FWHM. The results from the actual 3-yr WMAP temperature analysis are presented in a companion Letter.
Using a Gibbs sampling algorithm for joint CMB estimation and component separation, we compute the large-scale CMB and foreground posteriors of the 3-yr WMAP temperature data. Our parametric data model includes the cosmological CMB signal and instrum
ental noise, a single power law foreground component with free amplitude and spectral index for each pixel, a thermal dust template with a single free overall amplitude, and free monopoles and dipoles at each frequency. This simple model yields a surprisingly good fit to the data over the full frequency range from 23 to 94 GHz. We obtain a new estimate of the CMB sky signal and power spectrum, and a new foreground model, including a measurement of the effective spectral index over the high-latitude sky. A particularly significant result is the detection of a common spurious offset in all frequency bands of ~ -13muK, as well as a dipole in the V-band data. Correcting for these is essential when determining the effective spectral index of the foregrounds. We find that our new foreground model is in good agreement with template-based model presented by the WMAP team, but not with their MEM reconstruction. We believe the latter may be at least partially compromised by the residual offsets and dipoles in the data. Fortunately, the CMB power spectrum is not significantly affected by these issues, as our new spectrum is in excellent agreement with that published by the WMAP team. The corresponding cosmological parameters are also virtually unchanged.
