Sampling-based inference of the primordial CMB and gravitational lensing

The search for primordial gravitational waves in the cosmic microwave background (CMB) will soon be limited by our ability to remove the lensing contamination to $B$-mode polarization. The often-used quadratic estimator for lensing is known to be suboptimal for surveys that are currently operating and will continue to become less and less efficient as instrumental noise decreases. While foregrounds can, in principle, be mitigated by observing in more frequency bands, progress in delensing hinges entirely on algorithmic advances. We demonstrate here a new inference method that solves this problem by sampling the exact Bayesian posterior of any desired cosmological parameters, of the gravitational lensing potential, and of the delensed CMB maps, given lensed temperature and polarization data. We validate the method using simulated CMB data with nonwhite noise and masking on up to $650{\mathrm{deg}}^2$ patches of sky. A unique strength of this approach is the ability to perform joint inference of cosmological parameters, which control both the primordial CMB and the lensing potential, which we demonstrate here for the first time by sampling both the tensor-to-scalar ratio, $r$, and the amplitude of the lensing potential, $A_{\varphi }$. The method allows us to perform the most precise check to-date of several important approximations underlying CMB-S4 $r$ forecasting, and we confirm these yield the correct expected uncertainty on $r$ to better than 10%.

Figure 1

Typical simulated data and mask choices used in this work. Specifically, these correspond to the configuration BIG (see Table 2) with a true value of $r=0.01$. Reconstructed maps from this exact data are shown in Fig. 6. We note that an apodized pixel mask and an isotropic Fourier mask are not algorithm requirements, rather arbitrary choices we made for this example.

Figure 2

Three figures which are helpful in understanding the benefit of the reparametrization (described in Sec. 3), which makes sampling possible. Data configuration 2PARAM (see Table 2) is assumed for these figures. Left two panels: The reparamaterization includes switching from sampling the unlensed CMB fields, $f$, to sampling the “mixed” fields, $f^{\prime }$. The left two panels show that the power variation in a typical $f^{\prime }$ unmixed by various $\varphi$ is very small, an indication that large moves are allowed in Gibbs samples from the conditional $\mathcal{P}({\varphi }^{\prime }|f^{\prime },d)$. For comparison, the much larger variation in a typical $\tilde{f}$ when delensed by various $\varphi$ is shown, indicating that the lensed parametrization performs very poorly for polarization as is the case here. Right panel: Empirically, one finds that mixed $E$ is mostly lensed $E$ at all scales, while mixed $B$ is lensed $B$ at large scales but unlensed $B$ at small scales. We demonstrate this here by cross-correlating the mixed with the lensed or unlensed fields. This qualitatively conforms to the expectations of what should give a parametrization that is minimally degenerate (see discussion in Sec. 3).

Figure 3

Samples of $r$ and $A_{\varphi }$ at each iteration for two chains with the same data and starting random seed, but different choices for parametrizing the posterior (see Sec. 3). The blue line corresponds to using our optimized $\mathbb{G}(A_{\varphi })$ reparametrization, whereas the orange line shows the highly suboptimal choice of $\mathbb{G}(A_{\varphi })=\mathbb{1}$. No burn in is removed in either case. The simulated data here are generated according to configuration 2PARAM (see Table 2). The run time for a chain of this length is 48 hours on one GPU.

Figure 4

Posterior distribution for $r$ and $A_{\varphi }$ from a chain on simulated data in configuration 2PARAM (see Table 2). The samples that comprise this plot are shown in Fig. 3. For demonstration, here we use the getdist [63] plotting package rather than our Blackwell-Rao posterior density estimate. The ability to examine joint constraints on these parameters while performing optimal delensing for a realistic data set with masking is a unique strength of our approach. Here, we find these two parameters are highly uncorrelated, providing evidence that $A_{\varphi }$ can be fixed without impacting $r$ estimation. The orange curve shows a suboptimal choice of the $\mathbb{G}$ matrix, which causes that chain to converge more slowly.

Figure 5

Top panel: The trace of $r$ samples from chains in configuration BIG (see Table 2). Three different fiducial values of $r$ are explored, with the true value given by the black dashed line and each chain in a different color. No burn in is removed. Bottom panel: The same samples binned into histograms, as well Blackwell-Rao estimates of the posterior density, as described in Sec. 4d. These estimates recover very smooth distributions, even in the case where the true $r$ is zero and the constraint is just an upper bound.

Figure 6

True input maps (top row) as compared to posterior mean maps (bottom row) computed by averaging over chain samples. This chain uses configuration BIG (see Table 2) with a true value of $r=0.01$. The data for this chain are shown in Fig. 1. The first column shows the $\varphi$ map multiplied in Fourier space by $\ell$ to visually enhance smaller scales, the middle column shows $E$ modes that have been lensed into $B$, and the final column shows the reconstructed primordial $B$ modes.

Figure 7

Left column: Blackwell-Rao posteriors from each of ten chains on different simulated data for three different true values of the tensor-to-scalar ratio indicated in each row. The gray band is the product of the posteriors for each case, with the prior on $r$ importance sampled to be uniform (and with arbitrary normalization constant so as to fit on these axes). We expect that the gray band covers the fiducial value of $r$ to within its own width, as is indeed the case. This is a test of the coverage of our $\mathcal{P}(r|d)$ posteriors and hence, a test of the correctness of our procedure. Right column: The chain autocorrelation function for each of the chains in the left panel. The integrated autocorrelation time for these chains ranges from 5–33.

Figure 8

Left column: In blue, we overlay the power spectra of chain samples of $\varphi$ and of two quantities derived from $f$. The black dashed line gives the power spectrum of the truth, and the green line is the power spectrum of the posterior mean map. The three rows correspond to $\varphi$, unlensed $B$, and $E$-lensed-into-$B$. The posterior mean maps exhibit Wiener-filter like suppression, as expected, while the samples scatter around the true spectrum and quantify uncertainty. Right column: The same power spectra that are overlayed on the left but picking some specific multipoles and plotting the trace of their value throughout the chain. Visually one can see the acceptable correlation length of the power spectrum samples, as well as that they cluster around the input truth, confirmation that the likelihood is dominating over the prior that would otherwise pull these quantities towards zero. This is the same chain in configuration 2PARAM (see Table 2) used in Figs. 3 and 4.

Figure 9

Chain samples of the real part of the ten largest-scale and ten smallest-scale Fourier modes of the posterior $\varphi$, $B$, and $E$-lensed-into-$B$ maps. Each set of samples is normalized to unit variance, but the relative distance to the truth (shown in the black dashed line) is preserved. This chain uses configuration 2PARAM (see Table 2). Visually, we achieve great convergence even at the individual mode level. Out of the $\sim 200,000$ modes which are sampled, the only exceptions are perhaps the two largest scale $\varphi$ modes, which could benefit from a slightly longer chain. However, these two modes are not informative for $A_{\varphi }$, which remains very well converged (Fig. 3).

Figure 10

A comparison of different methods for forecasting constraints on $r$, assuming configuration MANY (see Table 2). Three different possible true values of $r$ are explored, indicated by vertical dotted lines. The dashed lines show expected Gaussian constraints forecasted with a method very similar to that used for CMB-S4. In solid lines, we show Gaussian distributions with standard deviation computed from the Fisher information on $r$. This work is the first time the Fisher information on $r$ from lensed CMB data has been calculated without approximation. In filled contours, we show the geometric mean of the posteriors from several chains. The excellent agreement between all of these is an important validation of the CMB-S4 $r$ forecasting methodology even in the presence of instrumental effects and masking, as is considered here.

Figure 11

The wall time in seconds for each Gibbs pass, for each configuration (see Table 2), and for running on CPU vs GPU. The CPU benchmarks utilize a full NERSC Cori Haswell node (Intel Xeon Processor E5-2698 v3), and the GPU benchmarks a single NVIDIA GTX 1080Ti GPU. Although our CPU code is highly optimized, our GPU code likely has room for significant improvement, despite already being faster.