Matrix based E/B separation for DC4 maps.

(updated 2019-02-17)   B. Racine, J. Willmert

In this posting, we report a first implementation of the purification matrix method for the S4 data challenge analysis.
(Update Feb 17th: added colorbars for the covariance matrix plots, and statistics of the power spectra of lensed-\(\Lambda\)CDM and noise simulations)

\( \newcommand{\mat}[1]{\mathbf{#1}} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\braket}[2]{\left\langle#1\middle|#2\right\rangle} \newcommand{\bb}{\tilde b} \newcommand{\vT}{\vphantom{\top}} \)

While full sky linear polarization fields can easily (and uniquely) be decomposed into E and B modes, cut-sky observations with non-uniform coverage and/or mode filtering have aliasing that provoke E to B leakages. The resulting bandpowers then contain pure-E, pure-B, and ambiguous modes that come from a mix of the original E and B modes. We can estimate the bias at the bandpower level by generating many simulations with a known power and computing the mean of these simulations' bandpowers compared to the input. But enhanced sample variance due to the leaked E-modes remains. More advanced estimators have been developed which avoid this increase in variance, by isolating the pure-B (or E) modes from the ambiguous modes.
In the data challenges, our maps are on a partial sky, and our noise is not homogeneous, which is modeled by dividing the noise maps by sqrt(hitcount map), where a hitcount maps depends on the pointing of our telescopes and the size of the field of view. In order to analyse the maps, we currently use the pure-s2hat (or p-s2hat) library, an implementation of the Smith pure estimator for cross-spectra, as introduced in this posting. As we saw in this short posting, while pure-s2hat reduces the E-to-B leakage at intermediate \(\ell\) as expected, there is a spurious low-\(\ell\) power. I think this is expected, coming from the 1/ell^2 dependence of the mixing kernels, and it could be computed and taken care of exactly. This excess of low-\(\ell\) power should be taken care of by our bandpower window functions, which are computed by imputing maps with delta function power of a red spectrum in the pipeline, and computing the bandpower response.

In this posting, we describe another power spectrum method, using matrix-based E/B separation.


The BICEP/Keck analysis pipeline uses a method described in the BKVII paper.

The main idea is to build two purification matrices that will project the pure-E and pure-B modes out of an observed map. In BKVII, it is shown that these purification matrices can be constructed by solving the following regularized generalized eigenvalue equation: \begin{align} \left( \mat{\tilde{C}_B} + σ^2\mat I \right) \mat x_i &= λ_i \left( \mat{\tilde{C}_E} + σ^2\mat I \right) \mat x_i , \label{eqn:geneig} \end{align} where \((λ_i,x_i)\) are the eigenvalues and corresponding eigenvectors, \(\mat{\tilde{C}_{(E/B)}}\) are the observed covariance matrices, \(\mat I\) is the identity matrix, and \(σ^2\) is the regularization parameter. The observed covariance matrix is \begin{align} \mat{\tilde{C}} = \mat R \mat C \mat R{}^\top , \label{eqn:obscovmat} \end{align} where the \(\mat R\) is the observation matrix, and \(\mat C\) is the covariance matrix. The observation matrix is here much simpler than in the BK case (where it includes all sorts of effects from data reduction). Here it should simply including the beam smoothing\(^*\), simple filterings\(^*\) and the pixel inverse variance weighting. The pixel-pixel covariance matrix is given by: \begin{align} \mat{C_{E}} \equiv \sum_{\ell m} C_\ell^{EE} \mat Y_{\ell m}^E \mat Y_{\ell m}^{E\dagger} , \label{eqn:covmat} \end{align} and similarly for \(\mat{C_{B}}\). The \(C_\ell\) is somewhat arbitrary here, but would ideally have a dependence similar to the signal. We here use a red spectrum: \(C_\ell = 1/\ell^2\). The covariance matrix is computed using the method from Tegmark & Oliviera-Costa (2001) , Appendix A. We focus only on the polarization.
\({}^*\): The beam smoothing and simple filterings are actually modeled at \(\ell\) level here, since it is a much simpler diagonal multiplication. We instead include it in the \(C_\ell\) used to compute the covariance matrix.

Once we have built these matrices, we can solve for equation \(\ref{eqn:geneig}\). The smallest eigenvalues correspond to the pure-E eigenmodes \(\mat e_i\) which are orthogonal to the B modes, whereas the highest eigenvalues correspond to the pure-B eigenmodes \(\mat b_i\) which are orthogonal to the E modes. The eigenvalues close to 1 are ambiguous modes.
We can then build the purification matrices, the description of which is in this posting (for those who have access to the BK pages). I will try to summarize it in an updated version of the current posting.
These are then applied to input observed maps \(\mat{\tilde{m}} \): \begin{align} \mat{\tilde{m}_{\rm Pure-E}}&= \mat{\Pi_{E} \tilde{m}}, \\ \mat{\tilde{m}_{\rm Pure-B}}&= \mat{\Pi_{B} \tilde{m}}, \\ \label{eqn:maps} \end{align}

The covariance matrix has a dimension \(N_{\rm pix} \times N_{\rm pix}\), so to make the computation tractable we here work with maps at \(N_{\rm side}=128\).

In section 1, we show the observation matrix for the different surveys (DC04 a,b,c and d)
In section 2, we show the observed covariance matrix, focusing on DC04c.
In section 3, we have a look at the eigensystem from equation \(\ref{eqn:geneig}\), and some eigenmodes.
In section 4, we build and plot the purification matrices (pure-E and pure-B projectors).
In section 5, we show pure-E and pure-B maps, as well as their power spectrum.
In section 6, we apply a low-pass filtering to the input maps to reduce high-\(\ell\) aliasing.

1. Observation matrix.

The observing matrix in the case of the S4 pipeline includes the effect from non-uniform coverage, i.e. the hitcount weighting in pixel space. Since we are also downsampling the maps to \(N_{\rm side}=128\), we have to include this in the observation matrix.

Figure 1 shows the different hitcount masks.
As explained here, 04 is the reference mask, with peak value at unity. The other masks have been rescaled to have the same total number of hits (i.e. same sum).
Note that these maps, at \(N_{\rm side}=128\), have respectively 10799, 68616, 9728 and 4915 observed pixels.

Figure 1:Hitcount maps: 04 is the circular 3% mask, 04b is the nominal Chile mask 04c is the nominal Pole mask and 04d is the hitcount used in BK14.

Figure 2 shows the different components of the observing matrix:

We also show the final observing matrix which is the product of the two.

Figure 2: Observation matrix (Q and U). The hitcount observing matrix is a \(\rm 2 N_{pix} \times 2 N_{\rm pix}\) matrix, where \(N_{\rm pix}\) is the number of observed pixels in the \(N_{side}=512\) hitcount map. The downsampling matrix is a \(\rm 2 N_{\rm pix} \times 2 N^{'}_{\rm pix}\) matrix where \(\rm 2 N^{'}_{pix} = 2 N_{pix}/16\). In this figure, we zoom on the first \(1000 \times 100\) elements, which only correspond to the Q map.

2. Covariance matrix.

We compute the covariance matrix using a EE and BB red spectra as described in \(\ref{eqn:covmat}\).
The covariance matrix is a \(\rm 2 N_{pix} \times 2 N_{pix}\) matrix. Column number n, represents the covariance of pixel number n with the rest of the pixels (both for Q and U).
In figure 3, we show the healpix projection of the column in the middle of the matrix.
Computing the full covariance matrix takes a while, though it's a fully parallelizable task. Some statistics of the runs on the oddyssey cluster are shown here.

Colorbars added on February 17th.. We plot the covariance matrix as is, with a free colorbar, as well as the covariance matrix divided by it maximum, with a symmetric colorbar.

Figure 3: Covariance matrix. The field button corresponds to the input red spectrum used to compute the covariance matrix.

Note that as a first test concerning the fact that we have high ell residuals in the current version with downsampling, we tried computing a covariance matrix where the input red spectrum stops at \(\ell_{\rm max} = 3*128-1=383\) instead of the \(\ell_{\rm max} = 700\) in the original run. In some of the following plots, we present results for both cases.

3. Solving the eigensystem.

Once we have the Q/U covariance matrices, both for E only and B only signals, we can solve the eigensystem shown in equation \(\ref{eqn:geneig}\), using \(\sigma^2 = 0.001 \).
Figure 4 shows the eigensystem, i.e. the sorted eigenvalue distribution, as well as the corresponding eigenvectors.
Figure 5 shows the first, fifth and hundredth eigen modes on the sky, as well as the modes corresponding to the 3 cutoff considered later on.

Figure 4: Eigensystem of the "04c" Pole case. The top plot shows the sorted eigenvalues, as well as the cutoff used later one when generating the pure-E and pure-B projectors. The bottom plot shows the corresponding eigenvectors, each column being normalized to have standard deviation of 1. The columns have a length equal to twice the number of observed pixels (for the Q and U part, and the line have a length equal to the number of eigenvalues.
The dotted, dashed and full vertical lines show the cutoff of 0.98, 0.9 or 0.5.
Figure 5: Eigenvector map corresponding to the eigenmode 1, 5, 100, and to the cutoff 0.5, 0.9, 0.98.

4. Pure-E and pure-B projectors.

We then build the purification matrices, or pure-E/pure-B projectors. We show such projectors in figure 6.

Figure 6: Pure-E and pure-B projectors. These are matrices with dimension \(2 N_{\rm pix} \times 2 N_{\rm pix} \) where \(2 N_{\rm pix}\) are the number of observed pixels. It is applied to the input Q-U vector map.

5. Maps and power spectra.

To test our purification matrices, we use the first unlensed \(\Lambda CDM\) map from the 04 data challenge. These maps don't contain any noise, and by definition only E modes. We expect to have a no power in the output pure-B map, and the input pure-E map to be similar to the input map. Note that we apodize these maps, as we would do in the case of the real analysis, using the inverse variance hitcount map: \begin{align} \mat{\tilde{m}_{\rm Pure-E/B}}&= \mat{\Pi_{E/B} (\tilde{m}*H)/H}., \\ \label{eqn:maps_apod} \end{align} where \(\mat H\) is the apodization map used. We have such maps plotted in figure 7. As we can see, when the beams are tight, we have some high ell residuals in the map.

Figure 7: Q and U maps before and after the projection of the pure-E and pure-B modes. The input map is the unlensed \(\Lambda\)CDM signal, so there should not be any B modes.

We then take the power spectra of these Q/U maps, using s2hat, (without purification, so healpix anafast should be giving similar results).
Figure 8 show the EE and BB power spectra for the different maps (input \(N_{\rm side}=512\) maps), pure-E and pure-B maps, as well as the pure-s2hat analysis of the input map and the downsampled input map (at \(N_{\rm side}=128\)).
It seems like the matrix purification is not efficient at high ell. But as can be seen for instance in this figure, s2hat behaves even worse when applied on a downsampled map.
Note also in this figure that the large beam map has a much better purification.
These two effects point towards high-\(\ell\) aliasing.
The two positive things to note are that when doing a fair comparison on the same downsampled input map, the matrix analysis seems to yield better results in terms of purification: first the residual B modes are lower, and the low-\(\ell\) residual in EE and BB is absent.

Figure 8: EE and BB Power spectra of the pure-E (dashed red) and pure-B maps (solid red), using (non pure) s2hat, compared to the input \(\rm{N_{side}}=512\) map's s2hat (dashed blue and pure-s2hat (dot-dashed blue) spectra. We also show in black the pure-s2hat results on the downsampled input map.

6. Low-pass filtering of the maps in preprocessing.

As we saw in the previous figures, there seems to be high-\(\ell\) aliasing in the power spectra.
As a first step to avoid this, we low-pass filter the maps. Here I use the original \(a_{\ell m}\) used for the map analysed before and multiply by a beam that has a roll-off at high-\(\ell\). I use a tapered cosine roll-off, for which I can control the \(\ell_{\rm min}\) and \(\ell_{\rm max}\) where the roll-off starts and ends.
We use a roll-off between 373 and 393, and one between 256 and 276.
As we see in Figure 9, the roll-off around \(3*N_{\rm side}-1\) is not strong enough, but the one at \(2*N_{\rm side}\) seems to remove the aliasing.

Figure 9: EE and BB Power spectra of the pure-E (dashed red) and pure-B maps (solid red), using (non pure) s2hat, compared to the input \(\rm{N_{side}}=512\) map's s2hat (dashed blue and pure-s2hat (dot-dashed blue) spectra. We also show in black the pure-s2hat results on the downsampled input map.

7. Pipelined analysis of 300 lensed-\(\Lambda\)CDM and noise simulations.

Added on February 17th.
In this section, we compare the mean and std of 300 power spectra using 3 methods: pure-s2hat, s2hat, and the purification matrix method. We analyze the maps after downsampling (and without the low-pass filtering from previous section).

Comments about the figure:

Figure 10: mean of the EE and BB power spectra of 300 downgraded maps, as well as the standard deviation.

8. Conclusion and perspective.

In this posting, we show the first implementation of the matrix-based purification method on the S4 maps. It is mostly reusing tools developed for the BK analysis.

A few natural follow-ups will be studied: