Wavemoth -- Fast spherical harmonic transforms by butterfly matrix compression

We present Wavemoth, an experimental open source code for computing scalar spherical harmonic transforms (SHTs). Such transforms are ubiquitous in astronomical data analysis. Our code performs substantially better than existing publicly available codes due to improvements on two fronts. First, the computational core is made more efficient by using small amounts of precomputed data, as well as paying attention to CPU instruction pipelining and cache usage. Second, Wavemoth makes use of a fast and numerically stable algorithm based on compressing a set of linear operators in a precomputation step. The resulting SHT scales as O(L^2 (log L)^2) for the resolution range of practical interest, where L denotes the spherical harmonic truncation degree. For low and medium-range resolutions, Wavemoth tends to be twice as fast as libpsht, which is the current state of the art implementation for the HEALPix grid. At the resolution of the Planck experiment, L ~ 4000, Wavemoth is between three and six times faster than libpsht, depending on the computer architecture and the required precision. Due to the experimental nature of the project, only spherical harmonic synthesis is currently supported, although adding support or spherical harmonic analysis should be trivial.Comment: 13 pages, 6 figures, accepted by ApJ

A multi-level solver for Gaussian constrained CMB realizations

Updated to match published version (no major changes)Updated to match published version (no major changes)Updated to match published version (no major changes)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/

Planck 2018 results. IV. Diffuse component separation

We present full-sky maps of the cosmic microwave background (CMB) and polarized synchrotron and thermal dust emission, derived from the third set of Planck frequency maps. These products have significantly lower contamination from instrumental systematic effects than previous versions. The methodologies used to derive these maps follow closely those described in earlier papers, adopting four methods (Commander, NILC, SEVEM, and SMICA) to extract the CMB component, as well as three methods (Commander, GNILC, and SMICA) to extract astrophysical components. Our revised CMB temperature maps agree with corresponding products in the Planck 2015 delivery, whereas the polarization maps exhibit significantly lower large-scale power, reflecting the improved data processing described in companion papers; however, the noise properties of the resulting data products are complicated, and the best available end-to-end simulations exhibit relative biases with respect to the data at the few percent level. Using these maps, we are for the first time able to fit the spectral index of thermal dust independently over 3 degree regions. We derive a conservative estimate of the mean spectral index of polarized thermal dust emission of beta_d = 1.55 +/- 0.05, where the uncertainty marginalizes both over all known systematic uncertainties and different estimation techniques. For polarized synchrotron emission, we find a mean spectral index of beta_s = -3.1 +/- 0.1, consistent with previously reported measurements. We note that the current data processing does not allow for construction of unbiased single-bolometer maps, and this limits our ability to extract CO emission and correlated components. The foreground results for intensity derived in this paper therefore do not supersede corresponding Planck 2015 products. For polarization the new results supersede the corresponding 2015 products in all respects

A CMB Gibbs sampler for localized secondary anisotropies

As well as primary fluctuations, CMB temperature maps contain a wealth of additional information in the form of secondary anisotropies. Secondary effects that can be identified with individual objects, such as the thermal and kinetic Sunyaev-Zel'dovich (SZ) effects due to galaxy clusters, are difficult to unambiguously disentangle from foreground contamination and the primary CMB however. We develop a Bayesian formalism for rigorously characterising anisotropies that are localised on the sky, taking the TSZ and KSZ effects as an example. Using a Gibbs sampling scheme, we are able to efficiently sample from the joint posterior distribution for a multi-component model of the sky with many thousands of correlated physical parameters. The posterior can then be exactly marginalised to estimate properties of the secondary anisotropies, fully taking into account degeneracies with the other signals in the CMB map. We show that this method is computationally tractable using a simple implementation based on the existing Commander component separation code, and also discuss how other types of secondary anisotropy can be accommodated within our framework

Multi-resolution Bayesian CMB component separation through Wiener filtering with a pseudo-inverse preconditioner

We present a Bayesian model for multi-resolution component separation for cosmic microwave background (CMB) applications based on Wiener filtering and/or computation of constrained realizations, extending a previously developed framework. We also develop an efficient solver for the corresponding linear system for the associated signal amplitudes. The core of this new solver is an efficient preconditioner based on the pseudo-inverse of the coefficient matrix of the linear system. In the full sky coverage case, the method gives an increased speed of the preconditioner, and it is easier to implement in terms of practical computer code. In the case where a mask is applied and prior-driven constrained realization is sought within the mask, this is the first time full convergence has been achieved at the full resolution of the Planck data set