Optimal Bayesian delensing progress update
Marius Millea, Dec 9, 2018
I'm working on applying our Bayesian delensing pipeline (from https://arxiv.org/abs/1708.06753 and another soon-to-be-out paper) to produce high fidelity r forecasts with optimal (beyond quadratic estimator) delensing, and to plug into the existing r forecasting pipeline. This post gives a first update of where I am at and what we are working on in collaboration with Ethan Anderes and Ben Wandelt.
Summary: I compute the joint maximum a posteriori estimator (ϕ_J) of the 02 polarization simulations, projected in the flat sky. These reconstructions have better (lower) variance than the quadratic estimate, but they're still not as good as they could be due to the mean-field induced by the flat-sky projection. I conclude we have to work in the curved sky, and doing so will be the subject of a followup post.
1. Background info on optimal delensing
Optimal delensing can most easily be understood from the joint Bayesian posterior probability:
- Eq 1: P(f,ϕ,θ|d)
where the f are the (T,E,B) fields, ϕ the lensing potential, θ any cosmological parameters, and d the data.
Our goal is to sample from this distribution, which yields a set of samples ϕᵢ and fᵢ (and correspondingly would yield a set of a samples of the lensing template f̃ᵢ = L(ϕᵢ)*fᵢ which can be used in the BK delensing pipeline). These samples are "optimal" in the sense that they use all information available, unlike the quadratic estimator which only uses quadratic combinations of the data and assumes the gradient approximation for lensing.
If you prefer thinking about estimators, a general statistical fact is that the posterior mean is the "optimal" estimator (ie minimizes the mean-squared error). So these samples produce the "optimal" estimate because the average over all the samples *is* the posterior mean.
Other estimators can be constructed which may approximate well the posterior mean. In Julien and Antony's paper (https://arxiv.org/abs/1704.08230), they take the posterior P(f,ϕ,θ|d), marginalize over f analytically,
- Eq 2: P(ϕ,θ|d) = ∫ df P(f,ϕ,θ|d)
and then find the best-fit ϕ, at fixed θ. We call this the ϕ_M estimator (the "M" referring to the "marginal" posterior which is maximized). We can also maximize the joint posterior distribution (Eq 1) directly, producing f_J and ϕ_J estimators (the "J" here standing for joint). ϕ_M has the advantage that the mask mean-field is not present, although ϕ_J is computationally easier to compute, and can be mean-field subtracted post-facto the same as the quadratic estimate. Its likely that ϕ_M is a lower variance estimator than ϕ_J (would be nice to check this), but both are much better than the quadratic estimate at these noise levels. The rest of this post computes ϕ_J, in preparation for sampling the joint distribution (where it seems computationally intractable to sample Eq 2 vs. Eq 1 for reasons I'd be happy to explain).
2. Can we use the flat-sky approximation for the deep patch?
A first question I wanted to answer is whether we can get away with the flat-sky aproximation for the deep (3% sky) patch. Although the curved-sky analysis is certaintly doable, being able to do it in flat-sky would make the analysis alot faster / simpler. To answer this, I did the following:
- Take the 02 set of polarization simulations (cmbs4_02_llcdm_f145_b04_ellmin30_map_2048_mc_*.fits), and project onto the flat sky (Lambert equal-area projection, 3.5' pixels).
- Produce a set of similarly pixelized simulations but directly on the flat-sky (ie assuming different k-modes are uncorrelated).
- Compute the ϕ_J estimate on each simulation, under the flat-sky approximation (i.e. assuming the k modes are uncorrelated, which is only true for the flat-sky sims). I use the same assumptions about the noise and instrument as in Lensing_reconstructions_02.00.
- Compute the total mean-field (which will include both mask and sky curvature) by averaging over the simulations.
- Compare the correlation coefficient of the mean-field subtracted reconstructions with the truth for the curved and flat-sky simulations.
Fig 1: First some basic checks on the convergence of the iterative ϕ_J procedure. Left figure shows that the χ² for a typical curved-sky run is asymptoting which is good, although it does so far away from the Gaussian expectation (the gray band). An example flat-sky sim shows the asymptote is usually much closer. So the curved-sky run is converged, but to a noticably bad fit (ie the model can't fit the sky curvature exactly right). The right figure shows in colored bands the power spectrum of the update to the current ϕ; by the final iteration we are making only very small fractional changes at relevant scales.
Fig 2: Here's what the mean-field subtracted reconstructed ϕ's look like. Its pretty cool to see by eye the improvement over the quadratic estimate. Even though the sky curvature means the model isn't exactly right, the mean field subtraction still forces the estimator to be unbiased, as suggested in this figure, just not necessarily the best possible variance. Note, the maximization procedure naturally produces a "wiener filtered" ϕ_J, hence the comparison against the wiener filtered quadratic estimate. Also ℓ<10 is filtered out in all images below for visual clarity.
Fig 3: Finally, we plot the correlation against the truth. There's improvement over the quadratic estimate, as expected. However, if we compare to the correlation we see from a flat-sky simulation, it falls short, which suggests we could get a better estimator by correctly modeling the sky curvature.
Note: The reconstructed ϕs from this post exist at NERSC and I'm happy to point any interested parties to them, although I'm not sure they're of much use yet beyond the excercise performed here.
3. Next steps
This is a first step towards getting samples from P(f,ϕ,θ|d), which I view as a valuable (but the only) way to do the S4 delensing analysis. I was hoping we could simplify things by working in the flat-sky approximation, but this post suggests that's not the case. Nevertheless, we've been working on computing ϕ_J on the curved sky, and it should be a matter of a few weeks, with sampling a few more after that. At that point, I could hand off posterior samples of lensed B (i.e. "delensing templates") which could be used in the BK style analysis.