Wiener filtering and pure E/B decomposition of CMB maps with anisotropic correlated noise

TL;DR

This paper introduces an advanced algorithm for reconstructing pure E and B mode maps from CMB data, effectively handling complex noise models and enabling statistically optimal separation crucial for detecting primordial B-modes.

Contribution

The paper presents an improved dual messenger algorithm that accounts for realistic anisotropic correlated noise and includes a method to estimate noise covariance from simulations.

Findings

Successfully reconstructs pure E and B maps free from mode leakage

Demonstrates high-speed execution suitable for large CMB datasets

Enables optimal Bayesian analysis for upcoming high-resolution CMB experiments

Abstract

We present an augmented version of our dual messenger algorithm for spin field reconstruction on the sphere, while accounting for highly non-trivial and realistic noise models such as modulated correlated noise. We also describe an optimization method for the estimation of noise covariance from Monte Carlo simulations. Using simulated Planck polarized cosmic microwave background (CMB) maps as a showcase, we demonstrate the capabilities of the algorithm in reconstructing pure E and B maps, guaranteed to be free from ambiguous modes resulting from the leakage or coupling issue that plagues conventional methods of E/B separation. Due to its high speed execution, coupled with lenient memory requirements, the algorithm can be optimized in exact global Bayesian analyses of state-of-the-art CMB data for a statistically optimal separation of pure E and B modes. Our algorithm, therefore, has a potentially key role in the data analysis of high-resolution and high-sensitivity CMB data, especially with the range of upcoming CMB experiments tailored for the detection of the elusive primordial B-mode signal.

Figures

Equations

$$\mathbf{d}=\mathcal{R}\mathbf{s}+\mathbf{n},$$

$${\chi }^2=(\mathbf{d}-\mathcal{R}\mathbf{s}{)}^{†}{\mathbf{N}}^{-1}(\mathbf{d}-\mathcal{R}\mathbf{s})+{\mathbf{s}}^{†}{\mathbf{S}}^{-1}\mathbf{s},$$

$${\mathbf{s}}_{\text{WF}}$$

$$\begin{gathered} {\chi }_{\mathbf{T},\mathbf{U}}^2=(\mathbf{d}-\mathbf{t}{)}^{†}{\bar{\mathbf{N}}}^{-1}(\mathbf{d}-\mathbf{t})+(\mathbf{t}-\mathcal{Y}\mathbf{B}\mathbf{u}{)}^{†}{\mathbf{T}}^{-1}(\mathbf{t}-\mathcal{Y}\mathbf{B}\mathbf{u}) \\ +(\mathbf{u}-\mathbf{s}{)}^{†}{\mathbf{U}}^{-1}(\mathbf{u}-\mathbf{s})+{\mathbf{s}}^{†}{\bar{\mathbf{S}}}^{-1}\mathbf{s}, \end{gathered}$$

$$\mathbf{u}$$

$$\mathbf{t}$$

$$\mathcal{A}={\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U}){\left[\mathfrak{R}{\mathbf{B}}^{†}{\mathcal{Y}}^{†}\mathcal{Y}\mathbf{B}{\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U})+\alpha 1\right]}^{-1}\mathfrak{R},$$

$$\tilde{\mathcal{A}}={\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U}){\left[\beta \mathfrak{R}{\mathbf{B}}^{†}\mathbf{B}{\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U})+\alpha 1\right]}^{-1}\mathfrak{R},$$

$${\mathcal{A}}^{-1}={\mathfrak{R}}^{†}\left[\mathfrak{R}{\mathbf{B}}^{†}{\mathcal{Y}}^{†}\mathcal{Y}\mathbf{B}{\mathfrak{R}}^{†}+\alpha (\bar{\mathcal{S}}+\mathbf{U}{)}^{-1}\right]\mathfrak{R}.$$

$${\mathbf{s}}_{n+1}={\mathbf{s}}_n+\tilde{\mathcal{A}}(\mathbf{b}-{\mathcal{A}}^{-1}{\mathbf{s}}_n),$$

$${\mathbf{N}}_i=\left(\begin{array}{ccc}⟨\mathcal{II}⟩& ⟨\mathcal{I}\mathcal{Q}⟩& ⟨\mathcal{I}\mathcal{U}⟩\\ ⟨\mathcal{Q}\mathcal{I}⟩& ⟨\mathcal{QQ}⟩& ⟨\mathcal{Q}\mathcal{U}⟩\\ ⟨\mathcal{U}\mathcal{I}⟩& ⟨\mathcal{U}\mathcal{Q}⟩& ⟨\mathcal{UU}⟩\end{array}\right),$$

$$\bar{\mathbf{N}}$$

$$=\mathbf{\Sigma }{\mathbf{P}}^{†}(\mathbf{\Delta }-\alpha \mathbf{P}{\mathbf{\Sigma }}^{-2}{\mathbf{P}}^{†})\mathbf{P}\mathbf{\Sigma },$$

$${\bar{\mathbf{N}}}^{-1}={\mathbf{\Sigma }}^{-1}{\mathbf{P}}^{†}(\mathbf{\Delta }-\alpha \mathbf{P}{\mathbf{\Sigma }}^{-2}{\mathbf{P}}^{†}{)}^{-1}\mathbf{P}{\mathbf{\Sigma }}^{-1},$$

$${\mathbf{t}}_{\mathbf{x}}={\left({\bar{\mathbf{N}}}^{-1}+{\mathbf{T}}^{-1}\right)}_{\mathbf{x}}^{-1}\left({\mathbf{T}}_{\mathbf{x}}^{-1}\mathcal{Y}\mathbf{B}{\mathbf{u}}_{\mathbf{\ell }}+{\bar{\mathbf{N}}}_{\mathbf{x}}^{-1}{\mathbf{d}}_{\mathbf{x}}\right).$$

$${\mathbf{s}}_{\text{CR}}=\stackrel{\sim }{\text{mathpzc}W}(\mathbf{d}-\hat{\mathbf{d}})+\hat{\mathbf{s}},$$

$$\mathbf{N}=\mathcal{Y}\mathbf{F}{\mathcal{Y}}^{†}\mathbf{D}\mathcal{Y}\mathbf{F}{\mathcal{Y}}^{†},$$

$${\mathbf{s}}_{\text{WF}}=\mathbf{S}(\mathbf{S}+\mathbf{F}{\mathcal{Y}}^{†}\mathbf{D}\mathcal{Y}\mathbf{F}{)}^{-1}{\mathcal{Y}}^{-1}\mathbf{d}.$$

$$\tilde{\mathbf{d}}=(\mathcal{Y}\mathbf{F}{\mathcal{Y}}^{†}{)}^{-1}\mathbf{d}.$$

$${\tilde{\mathbf{s}}}_{\text{WF}}$$

$$=(\mathbf{F}{\mathcal{Y}}^{†}{)}^{-1}\mathbf{S}(\mathbf{S}+\mathbf{F}{\mathcal{Y}}^{†}\mathbf{D}\mathcal{Y}\mathbf{F}{)}^{-1}{\mathcal{Y}}^{-1}\mathbf{d}$$

$$=(\mathcal{Y}\mathbf{F}{\mathcal{Y}}^{†}{)}^{-1}{\mathbf{s}}_{\text{WF}},$$

$$\mathbf{N}=\mathbf{D}\mathcal{Y}\mathbf{C}{\mathcal{Y}}^{†}\mathbf{D},$$

$$C_{\ell }=\frac{{\sigma }_N^2}{\beta }\left[1+{\left(\frac{{\ell }_{\mathrm{knee}}}{\ell }\right)}^{{\alpha }_{\mathrm{knee}}}\right],$$

$$\mathbf{u}$$

$$\cdot {\mathbf{B}}^{†}{\mathcal{Y}}^{†}{\mathbf{D}}^{-1}(\mathcal{Y}\varphi {\mathcal{Y}}^{†}{)}^{-1}{\mathbf{D}}^{-1}\mathbf{t}$$

$$\mathbf{t}$$

$$\cdot \left[(\mathcal{Y}\varphi {\mathcal{Y}}^{†}{)}^{-1}{\mathbf{D}}^{-1}\mathcal{Y}\mathbf{B}\mathbf{u}+(\mathcal{Y}\mathbf{C}{\mathcal{Y}}^{†}-\mathcal{Y}\varphi {\mathcal{Y}}^{†}{)}^{-1}{\mathbf{D}}^{-1}\mathbf{d}\right],$$

$$\begin{gathered} \tilde{\mathbf{t}}\equiv {\mathbf{D}}^{-1}\mathbf{t}=\mathcal{Y}(\mathbf{C}-\varphi 1){\mathcal{Y}}^{†}(\mathcal{Y}\mathbf{C}{\mathcal{Y}}^{†}{)}^{-1}{\mathbf{D}}^{-1}\mathcal{Y}\mathbf{B}\mathbf{u} \\ +(\mathcal{Y}\varphi {\mathcal{Y}}^{†})(\mathcal{Y}\mathbf{C}{\mathcal{Y}}^{†}{)}^{-1}{\mathbf{D}}^{-1}\mathbf{d}, \end{gathered}$$

$$\mathbf{v}$$

$$\mathbf{u}$$

$$\begin{gathered} \mathbf{u}={\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U}){\left[\mathfrak{R}{\mathbf{B}}^{†}{\mathcal{Y}}^{†}{\mathcal{M}}^{-1}\mathcal{Y}\mathbf{B}{\mathfrak{R}}^{†}(\bar{\mathcal{S}}+\mathbf{U})+\varphi \omega 1\right]}^{-1} \\ \cdot \mathfrak{R}{\mathbf{B}}^{†}{\mathcal{Y}}^{†}{\mathcal{M}}^{-1}\mathbf{v}, \end{gathered}$$

$$\tilde{\mathcal{A}}=({\beta }^{-1}\mathcal{Y}){\mathbf{C}}^{-1}({\beta }^{-1}{\mathcal{Y}}^{†}),{\mathcal{A}}^{-1}=\mathcal{Y}\mathbf{C}{\mathcal{Y}}^{†}.$$

$$\mathcal{A}$$

Full text

††pubyear: 2019††pagerange: Wiener filtering and pure $\mathcal{E}/\mathcal{B}$ decomposition of CMB maps with anisotropic correlated noise–B

1 Introduction

Cosmological inference from observations of the cosmic microwave background (CMB) polarization necessitates the separation of the contributions of the gradient and curl (or $\mathcal{E}$ and $\mathcal{B}$) components of the polarization signal to the data. These scalar $\mathcal{E}$ and pseudo-scalar $\mathcal{B}$ modes correspond to the spin-2 analogues of curl-free and divergence-free vector fields, respectively, with a polarization map being represented as the sum of both components (Zaldarriaga & Seljak, 1997; Seljak & Zaldarriaga, 1997; Kamionkowski et al., 1997a, b). The next generation of CMB experiments is focused on measuring the polarization of the CMB, with the $\mathcal{E}$-mode power spectrum providing an independent probe of the scalar modes measured via the temperature anisotropies (e.g. Abazajian et al., 2016; Suzuki et al., 2016; Henning et al., 2018; Louis et al., 2017). The $\mathcal{B}$-mode component of the polarization is the focus of growing interest in the community. First, they are an independent confirmation of the lensing effect detected in the temperature and $\mathcal{E}$-mode anisotropies, as $\mathcal{B}$ modes are produced from the gravitational lensing of $\mathcal{E}$ modes by the dark matter distribution along the line of sight (Zaldarriaga & Seljak, 1998). These lensing-induced $\mathcal{B}$ modes have been observed by high-resolution ground-based CMB experiments (e.g. Hanson et al., 2013; POLARBEAR Collaboration et al., 2017). Second, and most importantly, the detection of larger angular scale $\mathcal{B}$ modes, directly sourced by primordial gravitational waves, remains a key but elusive objective of modern cosmology (e.g. Guzzetti et al., 2016).

Observations of CMB polarization have attracted much interest due to the significance of the cosmological information encoded (e.g. Hu & White, 1997; Hu & Dodelson, 2002; Hu, 2003). The inclusion of $\mathcal{E}$-mode polarization data in parameter inference pipelines allows us to derive more stringent constraints (e.g. Galli et al., 2014), whilst the scientific potential of $\mathcal{B}$-mode anisotropy observations is extremely promising. A measurement of $\mathcal{B}$-mode signal on large angular scales ($\ell\la 100$), after discarding the expected lensed signal, would be regarded as a direct validation of the inflationary paradigm as the precursor of this stochastic background of gravitational waves (Kamionkowski & Kovetz, 2016). The amplitude of this background would directly constrain the energy scale of inflation, thereby ruling out some inflationary models (e.g. Zaldarriaga et al., 1997; Kinney, 1998; Tegmark et al., 2000), while also constraining the reionization period (Zaldarriaga, 1997). This would dramatically improve our understanding of the very early Universe.

The decomposition of the $\mathcal{E}$ and $\mathcal{B}$ modes on a partial sky is highly non-trivial due to the induced leakage between the two modes. Masked regions produce a discontinuity at the edges of the map and this results in a mixing of $\mathcal{E}$ and $\mathcal{B}$ modes, yielding ambiguous modes. Such modes can be sourced by either $\mathcal{E}$ or $\mathcal{B}$ components and cannot be uniquely assigned (Zaldarriaga, 2001; Lewis et al., 2002; Lewis, 2003; Bunn, 2002a, b; Bunn et al., 2003; Bunn, 2011). This is a highly prevalent issue due to most observations being made on an incomplete sky or full-sky maps being subjected to additional masking to reduce foreground contamination. Since the $\mathcal{E}$-mode power spectrum is much larger than that of $\mathcal{B}$ modes, the ambiguous modes significantly increase the variance of the $\mathcal{B}$ modes, resulting in a spurious measurement of the $\mathcal{B}$-mode power spectrum. The detection of the inflationary gravitational waves is especially challenging due to their relatively small amplitude and hence, efficient methods for pure $\mathcal{E}/\mathcal{B}$ decomposition are essential for extracting the cosmological information from CMB polarization data.

Several approaches for pure $\mathcal{E}/\mathcal{B}$ decomposition are described in the literature. While some techniques yield real-space maps of the derivatives of the polarization maps (e.g. Kim & Naselsky, 2010; Zhao & Baskaran, 2010; Kim, 2011; Bowyer et al., 2011), others are limited to power spectrum estimation via the construction of an eigenbasis for the pure-ambiguous decomposition (e.g. Challinor & Chon, 2005; Smith, 2006a, b; Smith & Zaldarriaga, 2007; Grain et al., 2009; Alonso et al., 2019). Nevertheless, such approaches do not result in pure $\mathcal{E}$ and $\mathcal{B}$ maps in the real space and are computationally intensive. Ferté et al. (2013) provides a quantitative comparison of the efficiency of the above techniques for power spectrum reconstruction. Wavelet-based techniques have also been proposed, but they must be carefully adapted for the problem under investigation (Cao & Fang, 2009; Rogers et al., 2016; Leistedt et al., 2017).

Bunn & Wandelt (2017) (hereafter BW) have recently shown that the $\mathcal{E}/\mathcal{B}$ decomposition can be approached from a Wiener filtering (Wiener, 1949) viewpoint, resulting in faster implementation as compared to the above methods, while providing real-space maps of the $\mathcal{E}$ and $\mathcal{B}$ modes. Another key advantage of such an approach is that it can be naturally extended to treat more interesting cases such as providing $\mathcal{E}$ maps free from the temperature anisotropy contributions by accounting for temperature and polarization correlations.

In this work, we present an augmented version of our dual messenger algorithm (Kodi Ramanah et al., 2017, 2018) for pure $\mathcal{E}/\mathcal{B}$ decomposition on the sphere, based on the principle of the Wiener filter. We adapt the algorithm to encode the BW prescription for reconstruction of pure $\mathcal{E}/\mathcal{B}$ maps and naturally extend the dual messenger framework to account for complex and realistic noise models. We demonstrate the application of this enhanced algorithm, designated as dante (DuAl messeNger filTEr), on a simulated CMB polarization data set which emulates the features of the actual Planck data.

The remainder of this paper is structured as follows. In Section 2, we provide a brief description of the dual messenger algorithm and illustrate a Jacobi relaxation method to account for the non-orthogonality of spherical harmonic transforms. We describe how it can be augmented to deal with non-trivial noise covariance in Section 3. We subsequently illustrate how the algorithm can encode the BW prescription for pure $\mathcal{E}/\mathcal{B}$ reconstruction, followed by an outline of the numerical implementation in Section 4. We then present a new optimization scheme for the estimation of noise covariance from Monte Carlo simulations and showcase the capabilities of dante in reconstructing pure $\mathcal{E}$ and $\mathcal{B}$ maps from simulated Planck data in Section 5. Finally, we summarize our main findings in Section 6. In Appendix A, we provide a brief description of spherical harmonic transforms, as employed in this work. Appendix B outlines the main steps in the derivation of the essential dual messenger equations to account for anisotropic correlated noise.

2 Dual Messenger Algorithm

We briefly review the underlying framework of the dual messenger algorithm for Wiener filtering polarized CMB maps. A complementary description is provided in Kodi Ramanah et al. (2018).

2.1 Wiener filter

In linear data analysis, we often encounter the computation of the so-called Wiener filter on large and complex data sets. The Wiener filter originates from the following statistical problem. We assume our observed data set $\bm{d}$ to be a linear combination of the signal $\bm{s}$ with covariance $\mathbf{S}$ and noise $\bm{n}$ with covariance $\mathbf{N}$, as follows:

[TABLE]

where the signal and noise covariances are given by $\mathbf{S}=\langle\bm{s}\bm{s}^{\dagger}\rangle$ and $\mathbf{N}=\langle\bm{n}\bm{n}^{\dagger}\rangle$, respectively. $\bm{\mathcal{R}}$ is the complete response operator, with the inclusion of harmonic transforms, beam and mask effects, that models the instrument response to incoming signal. It effectively corresponds to the overall model of how the instrument converts, on average, an incoming signal $\bm{s}$ to the observed data $\bm{d}$, with the residual being the noise $\bm{n}$, while encoding the relevant physics.

The Wiener filter solution, $\bm{s}_{\text{\tiny{{WF}}}}$, is the optimal linear filter when the signal and noise are both Gaussian random fields. For a particular realization of the data, $\bm{s}_{\text{\tiny{{WF}}}}$ therefore maximizes the posterior probability distribution $\propto\exp(-\chi^{2}/2)$, or equivalently minimizes:

[TABLE]

leading to the Wiener filter equation,

[TABLE]

$\bm{s}_{\text{\tiny{{WF}}}}$ is the least-square optimal solution: the Wiener filter minimizes the mean-square deviations $\langle\bm{\varepsilon}^{\dagger}\bm{\varepsilon}\rangle$ of the reconstruction errors $\bm{\varepsilon}=\bm{s}_{\text{\tiny{{WF}}}}-\bm{s}$, averaged over all signal and noise realizations.

The Wiener filter (Wiener, 1949) is a particularly important signal reconstruction technique, with ubiquitous applications in cosmology and astrophysics (e.g. Elsner & Wandelt, 2013, and references therein). Computing the Wiener filter solution for large and complex data sets from modern experiments, however, is numerically challenging. Indeed, the first matrix inversion above is dense in all bases and lives in a high-dimensional space. This space has typically the size of the number of elements in $\bm{d}$, which for Planck maps is $\mathcal{O}(10^{9})$, when accounting for polarization components and the nine frequency bands. Due to the size of the covariance matrices scaling as the square of the number of data samples, the storage and processing of dense systems become numerically intractable. Traditional approaches of computing the Wiener filter rely on costly and highly non-trivial numerical schemes, such as preconditioned conjugate gradient (PCG) methods (Eriksen et al., 2004; Wandelt et al., 2004; Smith et al., 2007; Seljebotn et al., 2014, 2019; Puglisi et al., 2018; Papež et al., 2018), requiring a preconditioner to approximate the matrix inversion involved. These complex techniques suffer from various numerical limitations, as discussed extensively in Elsner & Wandelt (2013) and Kodi Ramanah et al. (2017), when dealing with high-dimensional data sets. Such PCG methods do have some merits, however, as they are conducive to fast convergence, provided an adequate preconditioner tailored to the specific problem is available. A recent work by Horowitz et al. (2018) was based on recasting the Wiener filtering problem as an optimization scheme and provides an alternative promising approach for dealing with complex noise models. Münchmeyer & Smith (2019) recently proposed another interesting approach where a neural network was trained to Wiener filter CMB maps.

We recently presented the dual messenger algorithm, an enhanced variant of the standard messenger technique developed by Elsner & Wandelt (2013), as a general purpose Wiener filtering tool, which surpasses its predecessor in execution time (Kodi Ramanah et al., 2017). We have demonstrated the efficiency and unconditional stability of the dual messenger technique in Wiener filtering high-resolution polarized CMB data with correlated noise (Kodi Ramanah et al., 2018) (hereafter KLW18). As a preconditioner-free approach, this method circumvents the limitations of conventional PCG solvers in dealing with inherently ill-conditioned systems encountered in typical CMB polarization problems, as illustrated in KLW18.

The messenger method has been implemented within an efficient Gibbs sampling framework for Bayesian large-scale structure inference (Jasche & Lavaux, 2015), with this Gibbs-messenger sampling scheme subsequently adapted for CMB gravitational lensing (Anderes et al., 2015) and cosmic shear analyses (Alsing et al., 2016; Alsing et al., 2017). The messenger technique is being further developed for data analysis involving dense noise covariance matrices (Huffenberger, 2018) and has emerged as a promising CMB map-making tool (Huffenberger & Næss, 2018). The dual messenger algorithm has also been implemented in the field of optical and information engineering. For instance, it has been adapted for the removal of atmospheric haze from images (Fu et al., 2018), demonstrating the versatility of the tool developed. This class of messenger methods can therefore be tailored to solve a range of Wiener filtering problems and is not limited to astrophysical and cosmological applications.

2.2 Dual messenger field on the sphere

Conceptually, the essence of the messenger methods lies in the introduction of an auxiliary field to mediate between the different bases where the signal and noise covariances, $\mathbf{S}$ and $\mathbf{N}$, can be most conveniently expressed as sparse matrices. The addition of this messenger field allows the Wiener filter equation to be rewritten as a set of algebraic equations that must be solved iteratively, circumventing the requirement of matrix inversions or preconditioners.

With respect to the standard messenger technique, where a messenger field $\bm{t}$ is introduced at the level of the noise, the dual messenger algorithm incorporates an extra messenger field $\bm{u}$, at the level of the signal, with corresponding covariances $\mathbf{T}$ and $\mathbf{U}$, resulting in the following augmented $\chi^{2}$:

[TABLE]

where $\bar{\mathbf{N}}=\mathbf{N}-\mathbf{T}$ with $\mathbf{T}=\alpha\mathbb{1}$, where $\alpha\equiv\mathrm{min}(\mathrm{diag}(\mathbf{N}))$, and $\bar{\mathbf{S}}=\mathbf{S}-\mathbf{U}$ with $\mathbf{U}=\nu\mathbb{1}$, where $\nu\equiv\mathrm{min}(\mathrm{diag}(\mathbf{S}))$. When dealing with polarization fields, $\mu$ and $\nu$ are actually $3\times 3$ matrices, corresponding to the temperature, $\mathcal{E}$ and $\mathcal{B}$ components. $\bm{\mathcal{Y}}$ and $\bm{\mathcal{Y}}^{\dagger}$ correspond to the basis operators (synthesis and analysis operators, respectively) for the spherical harmonic transforms, as described in Appendix A, while $\mathbf{B}$ indicates convolution with an instrument beam. In terms of physical significance, $\bm{t}$ corresponds to a homogeneous portion of the noise covariance while its counterpart $\bm{u}$ is associated with the signal covariance. Optimizing $\chi^{2}_{\bm{T},\bm{U}}$ yields the following system of equations to be solved iteratively:

[TABLE]

Note that this is the reduced system of equations, with one of the messenger fields made implicit. The general system of equations is described in more depth in KLW18. To improve convergence, we implement a similar scheme as in KLW18. We artificially truncate the signal covariance $\mathbf{S}$ to some lower initial value of $\ell_{\mathrm{iter}}$, corresponding to a covariance $\mu$. We subsequently vary $\mathbf{U}$ to bring $\mu\rightarrow\nu$, such that in the limit $\mu=\nu$, we have $\bm{u}=\bm{s}$ and we recover the usual Wiener filter equation (3) from the above system of equations (5) and (6). This is formally valid as long as $\mu$ and $\nu$ are block matrices over harmonic space. We may therefore exploit this degree of freedom to solve the temperature and polarization signals at different rates. The above cooling scheme leads to a redefinition of $\bar{\mathbf{S}}$ using the Heaviside function as $\bar{\mathcal{S}}=\Theta(\mathcal{S}-\mathbf{U})$, where $\mathcal{S}$ corresponds to the eigenvalues of $\mathbf{S}$, i.e. $\mathbf{S}=\mathfrak{R}^{\dagger}\mathcal{S}\mathfrak{R}$ and $\bar{\mathbf{S}}=\mathfrak{R}^{\dagger}\bar{\mathcal{S}}\mathfrak{R}$.

To implement such a hierarchical scheme, we vary $\mathbf{U}$ via a cooling scheme for $\bm{\xi}$, where $\bm{\xi}=\mathbf{B}^{\dagger}\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{Y}}\mathbf{B}\mathbf{U}+\mathbf{T}$. To obtain the appropriate Wiener filter solution, we need to reduce $\mu\rightarrow\nu=0$, due to the continuous mode of the signal, i.e. the zero eigenvalue of $\mathbf{S}$. The cooling scheme for $\bm{\xi}$ entails gradually reducing $\bm{\xi}$ by a constant factor and iterating until $\bm{\xi}\rightarrow\mathbf{T}$, thereby finally attaining $\mu=0$, as required. A quantitative description of the rationale underlying the above cooling scheme is presented in our previous work (Kodi Ramanah et al., 2017).

dante is implemented in Python and it makes use of the HEALPix111http://healpix.jpl.nasa.gov (Górski et al., 2005) library, in particular the Python wrapper healpy, to perform the spherical harmonic transforms (SHTs). HEALPix employs an equal area projection scheme, where the SHTs are quasi-orthogonal, i.e. $\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{Y}}\mathbb{1}\approx(N_{\mathrm{pix}}/4\pi)\mathbb{1}\equiv\beta\mathbb{1}$, where $N_{\mathrm{pix}}$ denotes the number of pixels. We account for this non-orthogonality of SHTs via efficient Jacobi relaxation schemes, as described in future sections. If an equidistant cylindrical projection on a grid is adopted for the SHTs (e.g. Muciaccia et al., 1997; Huffenberger & Wandelt, 2010; McEwen & Wiaux, 2011), the equations presented in this work are significantly simplified, as a result of $\beta=1$. dante employs the Numba222https://numba.pydata.org (Lam et al., 2015) compiler for Python arrays and numerical functions to yield high-performance functions for all the required matricial manipulations to boost execution speed. Numba generates optimized native code using the LLVM compiler (Lattner & Adve, 2004) infrastructure and is used to parallelize the array operations.

2.3 Non-orthogonality of spherical harmonic transforms

Unlike in the case of discrete Fourier transforms, the spherical harmonic synthesis and analysis operators, i.e. $\bm{\mathcal{Y}}$ and $\bm{\mathcal{Y}}^{\dagger}$, respectively, are not orthogonal and differ by more than a transposition and a scale factor. The quality of the approximation, $\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{Y}}\mathbb{1}\approx(N_{\mathrm{pix}}/4\pi)\mathbb{1}\equiv\beta\mathbb{1}$, depends on the $\ell_{\mathrm{max}}$, $N_{\mathrm{pix}}$ and spherical grid considered.

To account for the non-orthogonality of spherical harmonic transforms, we incorporate an internal Jacobi relaxation method (Jacobi, 1845; Saad, 2003) in dante to refine the solution in harmonic space (cf. equation (5)). To obtain $\mathbf{U}$, we need the eigenvalues of $\mathbf{S}$, i.e. $\mathbf{U}=\mathrm{min}(\mathcal{S})\mathbb{1}$, where $\mathfrak{R}^{\dagger}\mathcal{S}\mathfrak{R}=\mathbf{S}$. Equation (5) can be formulated as $\bm{s}=\bm{\mathcal{A}}\bm{b}$, where $\bm{\mathcal{A}}$ is given by:

[TABLE]

after including the basis transformations, and $\bm{b}=\mathbf{B}^{\dagger}\bm{\mathcal{Y}}^{\dagger}\bm{t}$.

An approximation to $\bm{\mathcal{A}}$ can be obtained as follows:

[TABLE]

after using the approximate orthogonality relation $\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{Y}}\mathbb{1}\approx\beta\mathbb{1}$. The application of the operator $\bm{\mathcal{A}}$ is not well-defined but we nevertheless can apply its inverse $\bm{\mathcal{A}}^{-1}$ to a vector, by applying the relevant operators sequentially,

[TABLE]

We therefore make use of $\bm{\mathcal{A}}^{-1}$ and $\widetilde{\bm{\mathcal{A}}}$, to obtain the solution for $\bm{s}$ via the following Jacobi iterations:

[TABLE]

where $n$ denotes the number of Jacobi iterations.

The term $(\bar{\mathcal{S}}+\mathbf{U})^{-1}$ poses a numerical predicament for the final truncation in the signal covariance, where $\bar{\mathcal{S}}=\mathcal{S}-\mathbf{U}$ and $\mathbf{U}=\bm{0}$, and we subsequently require the inversion of $\mathcal{S}$. We circumvent the corner case due to the zero eigenvalues of the continuous modes in $\mathcal{S}$ by imposing the following constraint on the subspace $\mathcal{V}$ where $\mathcal{S}=\bm{0}$: $(\bar{\mathcal{S}}+\mathbf{U})^{-1}|_{\mathcal{V}}=\mathcal{S}^{+}$. $\mathcal{S}^{+}$ is the pseudo-inverse, i.e. $\mathcal{S}\mathcal{S}^{+}=\Pi$, where $\Pi$ is a projector.

2.4 Incomplete sky coverage

CMB data analysis inevitably requires the treatment of masks, with many practical applications requiring that parts of the sky be masked out. For full-sky observations, this is mainly to avoid contamination from the galactic foreground emissions, thereby preventing spurious power spectra measurements. In the case of ground-based or sub-orbital CMB experiments with partial sky coverage, missing data are accounted for using masks.

We provide an outline of the general procedure for solving the messenger equation (6) when dealing with temperature and polarization masks. Here, we assume correlated noise, such that the noise covariance $\mathbf{N}$ has the following block-diagonal form for every pixel $i$:

[TABLE]

where $\mathcal{I}$, $\mathcal{Q}$ and $\mathcal{U}$ are the Stokes parameters. More complex noise models will be described in Section 3.

We compute the covariance $\mathbf{T}$ of the messenger field $\bm{t}$ as follows: $\mathbf{T}=\mathrm{min}(\mathcal{N})\mathbb{1}=\alpha\mathbb{1}$, where $\mathfrak{D}^{\dagger}\mathcal{N}\mathfrak{D}=\mathbf{N}$. The noise covariance $\mathbf{N}$ can be written as $\mathbf{N}=\bm{\Sigma}\mathbf{C}\bm{\Sigma}=\bm{\Sigma}\mathbf{P}^{\dagger}\bm{\Delta}\mathbf{P}\bm{\Sigma}$, following the orthonormal decomposition of $\mathbf{C}$, where $\bm{\Sigma}$ is a diagonal matrix with the eigenvalues $\{\sigma_{\mathcal{I}}^{i},\sigma_{\mathcal{Q}}^{i},\sigma_{\mathcal{U}}^{i}\}$ corresponding to the noise amplitudes for the $i$th pixel, with the orthonormal decomposition of $\mathbf{C}$ resulting in the diagonal matrix $\bm{\Delta}$. We then obtain $\bar{\mathbf{N}}=\mathbf{N}-\mathbf{T}$ as follows:

[TABLE]

where $\alpha$, as stated above, is the smallest eigenvalue of $\mathbf{N}$. To solve the messenger equation (6), we require the inverse $\bar{\mathbf{N}}^{-1}$,

[TABLE]

such that $\bar{\mathbf{N}}^{-1}$ has a block-diagonal structure in pixel space. We obtain the solution for the messenger field by simply evaluating equation (6) in pixel space,

[TABLE]

We implement the temperature and polarization masks by increasing the noise variance to infinity for masked pixels, or numerically by setting the inverse noise covariance to zero. This is achieved by setting $\bm{\Sigma}^{-1}|_{\mathrm{mask}}=\bm{0}$, subsequently ensuring that data from masked regions do not contaminate the messenger field.

2.5 Constrained realizations

For full-sky coverage with parts of the sky masked out, we still seek the Wiener filter solution under the mask, constrained by the observations on the edges of the mask and determined by the prior inside the masked region. While this proposed reconstruction does not correspond to the true solution, it has the correct statistical properties, i.e. correct signal covariance. The generation of such constrained realizations is relevant for many practical CMB applications, such as exact likelihood evaluations via Gibbs sampling. A complementary conceptual discussion of the rationale underlying constrained realizations is provided in KLW18.

To draw Gaussian constrained realizations of the CMB sky, we need to simulate a reference signal $\hat{\bm{s}}$ in accordance with the prior signal covariance assumed for the Wiener filter. We also require a simulated data set $\hat{\bm{d}}$ whose signal and noise properties correspond to that of the data model. We adapt dante to generate constrained realizations (e.g. using the algorithm of Hoffman & Ribak, 1991) in only one application of the Wiener filter via:

[TABLE]

where $\stackrel{{\scriptstyle\sim}}{{\smash{\mathpzc{W}}\rule{0.0pt}{4.73611pt}}}$ indicates the application of the dual messenger operator at a given precision and cooling scheme. We therefore only need to modify the input data fed to dante to draw unbiased constrained realizations of the signal that are consistent with the observed data, i.e. having the correct covariance properties.

3 Dual Messenger Generalizations

The dual messenger approach can be extended to a broader class of problems, accounting for highly non-trivial noise covariance, as described in the following sections. In practice, the structure of the noise covariance, in pixel space, is influenced by the noise properties of the CMB time-ordered data and the scanning pattern of the telescope.

3.1 Correlated modulated noise

The first possible generalization for the noise model is the case of correlated modulated noise. The noise covariance, in pixel space, can be written as:

[TABLE]

where $\mathbf{F}$ is a smoothing kernel which is diagonal in the same basis as $\mathbf{S}$, i.e. harmonic space, while $\mathbf{D}$ is the noise variance that can be easily diagonalized in some other basis, e.g. pixel space. The desired Wiener filter from equation (3) is then

[TABLE]

It turns out that this can be solved directly by the dual messenger scheme described above without any modifications. We transform the data via a simple pre-whitening step, as follows:

[TABLE]

The Wiener filter for this model can then be computed via the following steps:

[TABLE]

after plugging in the effective data vector given by equation (18).

This problem therefore reduces to one that can be solved directly using the dual messenger algorithm, requiring the same computational time as in the white noise case. The only additional steps required are a simple pre-whitening, followed by a post-smoothing operation with $\mathbf{F}^{-1}$ and $\mathbf{F}$, respectively. For this particular case, we do not demonstrate the application of dante, as the implementation is straightforward.

3.2 Modulated correlated noise

The second generalization of the noise model corresponds to modulating the amplitude of spatially correlated noise. This is a more realistic noise model, typical of CMB experiments such as Planck, resulting from the scanning strategy of the instrument. The noise covariance, in pixel space, now takes the following form:

[TABLE]

where $\mathbf{C}$ is the isotropic noise covariance, encoding the inverse frequency ($1/f$) noise correlation on the large scales, typically associated with atmospheric noise, and therefore diagonal in harmonic space. The modulation, described by $\mathbf{D}$, is sparse in pixel space. The power spectrum $C_{\ell}$ of the non-modulated part of the noise can be expressed as:

[TABLE]

with the characteristic scale $\ell_{\mathrm{knee}}$ and tilt $\alpha_{\mathrm{knee}}$ of the power-law component, with $\sigma_{N}$ being the noise amplitude per pixel.

The modulation takes into account the variation of noise amplitude due to the amount of integration time spent in any single pixel. The isotropic noise covariance in the map is expected to be a good approximation for scan strategies that cross each pixel in many directions, thereby isotropizing the way the time-ordered data noise correlations project onto the sky. But even in cases where the distribution of scan directions per pixel is not entirely isotropic, such as for the Planck data, this noise model was found to be of sufficient quality to derive the optimal primordial non-Gaussianity estimators (Planck Collaboration et al., 2016c, 2019).

As a result of the above dense noise covariance, equation (6) is no longer algebraically solvable due to the required inversion of a fully dense system. But since we are free to choose the covariance $\mathbf{T}$ of the messenger field $\bm{t}$, we set $\mathbf{T}=\mathbf{D}(\bm{\mathcal{Y}}\phi\bm{\mathcal{Y}}^{\dagger})\mathbf{D}$, where $\phi=\mathrm{min}(\mathrm{diag}(\mathbf{C}))$. This yields the following system of equations:

[TABLE]

where the second equation can be simplified to the following form via straightforward linear algebraic manipulations:

[TABLE]

which can now be solved trivially to obtain the messenger field. The first equation (22), however, cannot be solved directly, but it can be conveniently expanded using an extra messenger field $\bm{v}$, with covariance $\mathbf{V}=\omega(\bm{\mathcal{Y}}\phi\bm{\mathcal{Y}}^{\dagger})\mathbb{1}$, where $\omega\equiv\mathrm{min}(\mathrm{diag}(\mathbf{D}^{2}))$, yielding the following two trivially solvable equations:

[TABLE]

where the coupling matrix $\bm{\mathcal{M}}$ is defined as $\bm{\mathcal{M}}\equiv\bm{\mathcal{Y}}\bm{\mathcal{Y}}^{\dagger}$. We therefore must solve the above system of three equations (24)-(26) when accounting for modulated correlated noise. Equation (26) can be written explicitly in terms of the relevant basis transformations as follows:

[TABLE]

where, as before, $\mathbf{S}=\mathfrak{R}^{\dagger}\mathcal{S}\mathfrak{R}$. An in-depth derivation of these equations is laid out in Appendix B. We encode the mask by doing the decomposition, $\mathbf{D}=\bm{\Sigma}\widetilde{\mathbf{D}}\bm{\Sigma}=\bm{\Sigma}\mathbf{P}^{\dagger}\bm{\Delta}\mathbf{P}\bm{\Sigma}$, as described in Section 2.4, and setting $\bm{\Sigma}^{-1}|_{\mathrm{mask}}={\bm{0}}$. This prescription corresponds exactly to dropping the contribution of observations that are considered masked out, which may be deduced from equation (4). We apply the cooling scheme to $\bm{\xi}=\mathbf{U}+\phi\omega\mathbb{1}$, as described in Section 5.3.

3.3 Nested Jacobi relaxation

As mentioned in Section 2.2, the above equations (24), (25) and (27) may be simplified significantly using the approximation $\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{Y}}\mathbb{1}\approx\beta\mathbb{1}$, thereby reducing the required number of SHTs. This approximation is not exact due to the coupling of the SHTs on the pixelized sky. In this work, we employ Jacobi relaxation to correct for the operations $\bm{\mathcal{M}}^{-1}\equiv(\bm{\mathcal{Y}}\bm{\mathcal{Y}}^{\dagger})^{-1}$ and $(\bm{\mathcal{Y}}\mathbf{C}\bm{\mathcal{Y}}^{\dagger})^{-1}$.

We follow a similar rationale and employ the same notation as in Section 2.3, with the Jacobi iteration given by equation (10), where $\bm{b}$ is any arbitrary vector in harmonic space. For the case of $\bm{\mathcal{A}}=(\bm{\mathcal{Y}}\mathbf{C}\bm{\mathcal{Y}}^{\dagger})^{-1}$, the respective operations are as follows:

[TABLE]

The correction for $\bm{\mathcal{M}}^{-1}$ is analogous to the above, with $\mathbf{C}$ set to identity matrix, and must be embedded within the less trivial Jacobi relaxation for equation (27). The resulting nested relaxation scheme requires the following operations:

[TABLE]

and the corresponding inverse of operator $\bm{\mathcal{A}}$ given by:

[TABLE]

with the basis operations included, and $\bm{b}=\mathbf{B}^{\dagger}\bm{\mathcal{Y}}^{\dagger}\bm{\mathcal{M}}^{-1}\bm{v}$.

4 Pure $\mathcal{E}/\mathcal{B}$ decomposition via Wiener filtering

In a finite patch of sky, the polarization field cannot be uniquely decomposed into pure $\mathcal{E}$ and pure $\mathcal{B}$ modes. Nevertheless, the polarization map can be uniquely decomposed into three distinct components, commonly referred to as the “pure $\mathcal{E}$”, “pure $\mathcal{B}$” and “ambiguous” modes (Lewis et al., 2002; Bunn et al., 2003). This new set of ambiguous modes receives contributions from both $\mathcal{E}$ and $\mathcal{B}$ modes. In such a framework, the signal vector space is divided into three orthogonal subspaces. The pure $\mathcal{B}$ modes exist on the vector subspace orthogonal to that of all $\mathcal{E}$ modes, and similarly for the pure $\mathcal{E}$ modes. The ambiguous component, however, lies in the subspace orthogonal to both pure $\mathcal{E}$ and $\mathcal{B}$ subspaces. This decomposition ensures that a reconstructed pure $\mathcal{B}$ map is not contaminated by $\mathcal{E}$ modes.

$\mathcal{E}/\mathcal{B}$ separation methods based on this pure-ambiguous decomposition originally involved the construction of an eigenbasis for the various orthonormal subspaces, but this is a tedious and numerically slow procedure. The $\mathcal{E}/\mathcal{B}$ decomposition is trivial for exact methods such as Gibbs sampling (Larson et al., 2007), which infer the posterior statistics of a full-sky signal conditional on the data. Gibbs sampling requires a complete sky sample, i.e. optimally filtered data augmented to cover the whole sky via constrained generalizations (see also KLW18). This is the basis of the motivation behind the Wiener filtering approach proposed by Bunn & Wandelt (2017).

We briefly review the rationale and the formalism behind this new method, and describe how it can be incorporated in dante. A more comprehensive description is provided in Bunn & Wandelt (2017).

4.1 Background and notation

Considering only polarization measurements, the data set can be described as a $2N_{\mathrm{pix}}$ dimensional vector of the Stokes parameters $\mathcal{Q}$ and $\mathcal{U}$, i.e. $\bm{d}=(\bm{d}_{\mathcal{Q}},\bm{d}_{\mathcal{U}})$, where $N_{\mathrm{pix}}$ corresponds to the dimension of the pixelized map of a given Stokes parameter. We account for the contribution from the temperature anisotropy, Stokes $\mathcal{I}$, in a future section. We assume a data model as given by equation (1), and Gaussian white noise, although the results presented below would still be relevant for more complex noise covariance.

The signal can be expressed as a spherical harmonic expansion,

[TABLE]

where $\hat{r}_{j}$ labels the pixel corresponding to measurement $j$, while the index $\mathpzc{z}\in\{\mathcal{Q},\mathcal{U}\}$ denotes the associated Stokes parameter. The functions $Y$ can be expressed in terms of spin-2 spherical harmonics:

[TABLE]

The signal can therefore be written as

[TABLE]

with the coefficients $a_{\ell m}^{\mathcal{E}}$ and $a_{\ell m}^{\mathcal{B}}$ encoded in the vectors $\bm{e}$ and $\bm{b}$, respectively. The matrices $\mathbf{Y}_{\mathpzc{Z}}$, for $\mathpzc{Z}\in\{\mathcal{E},\mathcal{B}\}$, consist of the spherical harmonics evaluated at the given pixel locations.

Under the assumption of data sourced by a statistically isotropic and random Gaussian process, the signal is uniquely described by the following covariance:

[TABLE]

For a full-sky data set, where $\bm{d}$ covers the whole sky, the matrices $\mathbf{Y}_{\mathcal{E}}$ and $\mathbf{Y}_{\mathcal{B}}$ span orthogonal spaces, and hence the $\mathcal{E}$-$\mathcal{B}$ coupling issue does not arise and the decomposition is straightforward. Incomplete sky coverage, however, results in the ambiguous modes that lie in both subspaces at the same time. The pure $\mathcal{B}$ space can therefore be described as the orthogonal complement of the space spanned by $\mathbf{Y}_{\mathcal{E}}$, with an analogous definition for the pure $\mathcal{E}$ space. The signal vector space can consequently be divided into three orthogonal subspaces, with the third ambiguous space being orthogonal to both of the pure subspaces. By projecting the data vector $\bm{d}$ on the pure $\mathcal{B}$ subspace, the $\mathcal{E}$-mode signal is mapped onto the null space of $\mathbf{Y}_{\mathcal{B}}$, ensuring no contamination of $\mathcal{E}$ modes in the resulting pure $\mathcal{B}$ map.

4.2 Purification framework

If the signal covariance $\mathbf{S}$ employed in the Wiener filter equation (3) contains covariances of both $\mathcal{E}$ and $\mathcal{B}$ signals, then the resulting Wiener-filtered map would have contributions from both the scalar and pseudo-scalar components, i.e. $\bm{s}_{\text{\tiny{{WF}}}}=\bm{s}_{\text{\tiny{{WF}}}}^{\mathcal{E}}+\bm{s}_{\text{\tiny{{WF}}}}^{\mathcal{B}}$. Bunn & Wandelt (2017) proposed the following approach to isolate them from each other: Conceptually, the rationale is to treat one component as a source of noise. We can obtain the Wiener-filtered $\mathcal{B}$ map, for instance, via the following replacements: $\mathbf{S}\rightarrow\mathbf{S}_{\mathcal{B}}$ and $\mathbf{N}\rightarrow\mathbf{S}_{\mathcal{E}}+\mathbf{N}$, such that the Wiener filter equation (3) becomes:

[TABLE]

with an analogous expression for $\bm{s}_{\text{\tiny{{WF}}}}^{\mathcal{E}}$. Recall that $\bm{s}_{\text{\tiny{{WF}}}}=\bm{s}_{\text{\tiny{{WF}}}}^{\mathcal{E}}+\bm{s}_{\text{\tiny{{WF}}}}^{\mathcal{B}}$.

However, the above Wiener-filtered maps are contaminated by the ambiguous modes. Due to our prior signal covariance having a much higher $\mathcal{E}$-mode power, the theory assigns the ambiguous modes with high signal-to-noise mostly to the $\mathcal{E}$ map. In order to ensure no cross-contamination, whereby a pure $\mathcal{B}$ map should have contributions strictly from $\mathcal{B}$ modes, Bunn & Wandelt (2017) suggested raising the signal covariance associated to the $\mathcal{E}$ component to infinity, and provided a proof that this gives the same result as doing a costly eigenmode decomposition and projecting out the ambiguous modes. We define

[TABLE]

as the signal covariance with the $\mathcal{E}$-mode power amplified by a factor of $\lambda$. Substituting $\mathbf{S}_{\mathcal{B}}+\mathbf{S}_{\mathcal{E}}\rightarrow\mathbf{S}(\lambda)$ in equation (37) yields

[TABLE]

such that in the limit $\lambda\rightarrow\infty$, only the pure $\mathcal{B}$ modes survive in the null space of $\mathbf{S}_{\mathcal{E}}$. A strictly analogous procedure holds for the pure $\mathcal{E}$ component.

4.3 Numerical implementation

To facilitate the numerical evaluation of the expressions above, it is more convenient to work with a full-sky data set, so that we can use fast transforms to move back and forth between the pixel and spherical harmonic spaces. Masked pixels are assigned infinite noise covariance, i.e. we set the inverse noise covariance to zero.

Assuming isotropic and Gaussian CMB anisotropies, the signal covariance $\mathbf{S}$ is diagonal in spherical harmonic basis,

[TABLE]

with the ordering convention of having the $\mathcal{E}$-mode component first. Hence, we have

[TABLE]

As a result, $\mathbf{S}(\lambda)^{-1}|_{\lambda\rightarrow\infty}\rightarrow\mathbf{S}_{\mathcal{B}}^{+}$, where the pseudo-inverse $\mathbf{S}_{\mathcal{B}}^{+}$ is the inverse of $\mathbf{S}_{\mathcal{B}}$ within the subspace spanned by $\mathbf{Y}_{\mathcal{B}}$ and is zero in the orthogonal subspace of $\mathcal{E}$ modes. We consequently obtain the operator that projects onto the pure $\mathcal{E}$ subspace as $\mathbf{S}_{\mathcal{B}}\mathbf{S}_{\mathcal{B}}^{+}=\mathbf{P}_{\mathcal{B}}$. Finally, we obtain the Wiener-filtered pure $\mathcal{B}$ map as

[TABLE]

The above formalism can be generalized to include the correlation between polarization and temperature anisotropies, where the data, $\bm{d}=(\bm{d}_{\mathcal{I}},\bm{d}_{\mathcal{Q}},\bm{d}_{\mathcal{U}})$, are pixelized maps of the Stokes parameters, $\mathcal{I}$, $\mathcal{Q}$ and $\mathcal{U}$, where $\mathcal{I}$ here corresponds to the temperature anisotropy. The signal covariance matrix now has a block-diagonal structure in spherical harmonic space, with a $3\times 3$ sub-matrix, for all multipole moments $\ell$, as follows:

[TABLE]

with the vanishing cross-spectra, $C_{\ell}^{TB}$ and $C_{\ell}^{EB}$, set to zero.

To find the pure $\mathcal{B}$ map, we proceed as before, i.e. $\mathbf{S}_{\mathcal{E}}\rightarrow\lambda\mathbf{S}_{\mathcal{E}}$ with $\lambda\rightarrow\infty$. In this limit, the temperature and polarization components decouple, yielding the following signal covariance:

[TABLE]

such that the above equation (43) for the Wiener-filtered pure $\mathcal{B}$ map still holds.

A similar reasoning results in the following equation for the Wiener-filtered pure $\mathcal{E}$ map:

[TABLE]

where $\mathbf{S}_{\mathcal{TE}}$ is the signal covariance corresponding to:

[TABLE]

$\mathbf{S}_{\mathcal{TE}}^{+}$ is then the pseudo-inverse associated to the subspace containing the temperature and $\mathcal{E}$ modes, i.e. all the $\mathcal{B}$ modes lie in the null space of both $\mathbf{S}_{\mathcal{TE}}$ and $\mathbf{S}_{\mathcal{TE}}^{+}$.

In this work, we encode this prescription in dante for optimal reconstruction of pure $\mathcal{E}$ and $\mathcal{B}$ maps via Wiener filtering. The numerical implementation for the pure $\mathcal{B}$ case entails the usual Wiener filtering procedure, i.e. solving equation (3), but assuming infinite covariance for the $\mathcal{E}$ component, followed by the application of the relevant projection operator $\mathbf{P}_{\mathcal{B}}$ to obtain the pure $\mathcal{B}$ map. An analogous procedure yields the pure $\mathcal{E}$ map.

The above formalism still holds in the presence of more complex noise models, such as the anisotropic correlated noise considered here. Since the signal covariance becomes fully diagonal when reconstructing the pure $\mathcal{E}$ or $\mathcal{B}$ map, we solve equation (26) itself since no basis transformations are required, as follows:

[TABLE]

There is, nevertheless, a caveat in the implementation of the above equation. The signal covariance $\mathbf{S}$, as given by equations (45) and (47), is actually the pseudo-inverse of the well-defined $\mathbf{S}(\lambda)^{-1}$ in the limit $\lambda\rightarrow\infty$. The trivial implementation is to set the relevant components of $\mathbf{S}_{\mathcal{TE}}$ and $\mathbf{S}_{\mathcal{B}}$ to a numerically large value. Alternatively, equation (48) can be expressed in the following numerically convenient form:

[TABLE]

which, for the final step of the cooling scheme, i.e. $\mathbf{U}=\bm{0}$, reduces to:

[TABLE]

We verified that both implementations resulted in identical solutions, within the limit of numerical errors.

5 Numerical experiments

In this section, we demonstrate the application of dante to an artificially generated but realistic CMB polarization data set. We present the procedure for the mock generation, followed by a description of the different steps in the data analysis pipeline.

5.1 Mock generation

To simulate joint temperature and polarization maps on the sphere, we make use of healpy to generate realizations of $a_{\bm{\ell}}^{\mathcal{T}}$, $a_{\bm{\ell}}^{\mathcal{E}}$ and $a_{\bm{\ell}}^{\mathcal{B}}$ signals with the correct covariance properties (cf. equation (44)), taking into account the correlation between CMB temperature anisotropy and polarization. We employed camb333http://camb.info (Lewis et al., 2000) to generate the input angular power spectra, $C^{TT}_{\ell}$, $C^{EE}_{\ell}$, $C^{BB}_{\ell}$ and $C^{TE}_{\ell}$, from which the corresponding CMB signals are drawn. We assume a standard $\Lambda$CDM cosmology with the set of cosmological parameters ($\Omega_{\mathrm{m}}=0.32$, $\Omega_{\Lambda}=0.69$, $\Omega_{\mathrm{b}}=0.05$, $h=0.67$, $\sigma_{8}=0.83$, $n_{\mathrm{s}}=0.97$) from Planck (Planck Collaboration et al., 2016b). We can then construct the input $\mathcal{Q}$ and $\mathcal{U}$ maps by transforming realizations of $\mathcal{E}$ and $\mathcal{B}$ signals (cf. equation (64) in Appendix A) on the sphere, with HEALPix resolution of $N_{\mathrm{side}}=128$ and $\ell_{\mathrm{max}}=128$, such that the total number of pixels is $N_{\mathrm{pix}}=12\times N_{\mathrm{side}}^{2}\approx 2\times 10^{5}$. The input Stokes parameters’ maps are subsequently contaminated with modulated correlated noise, as described by equation (51) below, with a white noise amplitude of $\sigma_{N}=40.0\>\mu\mathrm{K}$ per pixel, typical of high-sensitivity CMB experiments tailored for the detection of $\mathcal{B}$ modes, and the corresponding $1/f$ noise parameters of $\ell_{\mathrm{knee}}=10$ and $\alpha_{\mathrm{knee}}=1.5$ (cf. equation (21)). We employ the smica Planck temperature and polarization masks (Planck Collaboration et al., 2016a), corresponding to sky fractions of $f_{\mathrm{sky}}^{\mathcal{T}}=0.76$ and $f_{\mathrm{sky}}^{\mathcal{P}}=0.84$, respectively, as depicted in Fig. 1. While our formalism and code account for the effect of a beam, we set the beam operator to identity for our present test cases.

5.2 Estimation of noise covariance

We now present a posterior optimization method to estimate the noise covariance using Monte Carlo (MC) simulations. For the case of modulated correlated noise, the data can be modelled as follows:

[TABLE]

following the notation of equation (1), where $\mathbf{D}$ is diagonal in pixel space. We also use liberally the notation $\mathbf{C}^{1/2}$ to indicate the positive square root matrix of $\mathbf{C}$. As described in Section 3.2, the noise covariance is now given by

[TABLE]

with $\mathbf{C}$ being the isotropic, homogeneous, noise covariance, which incorporates the inverse frequency ($1/f$) noise correlation on the large scales. The overall aim is to estimate $\mathbf{D}$ and $\mathbf{C}$ using MC simulations by casting the covariance estimation problem as a two-level optimization scheme. The noise realizations can be modelled as the following linear combination:

[TABLE]

where $\bm{z}$ and $\bm{k}$ are Gaussian random fields with covariances, $\langle\bm{z}\bm{z}^{\dagger}\rangle=\mathbb{1}$ and $\langle\bm{k}\bm{k}^{\dagger}\rangle=\kappa^{2}\mathbb{1}$, respectively. The corresponding $\chi^{2}$, as the negative of the logarithm of the posterior distribution, with the sum over the contribution of each MC realization, can be written as:

[TABLE]

To obtain the maximum a posteriori estimate of $\mathbf{D}$ and $\mathbf{C}$, we must optimize the above $\chi^{2}$ with respect to $\bm{z}_{i}$ and $\mathbf{D}$, in the limit $\kappa\to 0$. The $\chi^{2}$ optimization with respect to $\bm{z}_{i}$ yields, for a given MC simulation,

[TABLE]

which, in the limit $\kappa\to 0$, simplifies to

[TABLE]

where $\tilde{\Pi}^{\dagger}=(\bm{\mathcal{Y}}\bm{\mathcal{Y}}^{-1})^{\dagger}$ is a projector onto the pixel subspace. The inversion per matrix is acceptable for the term in parenthesis because the operation $\bm{\mathcal{Y}}^{\dagger}$ already projects on the subspace of maps bandwidth limited to $\ell_{\text{max}}$. Within that space the $\bm{\mathcal{Y}}$ operator becomes invertible, though at some cost. Optimizing the $\chi^{2}$ with respect to $\mathbf{D}$ leads to

[TABLE]

where we defined $\bm{\mathcal{E}}\equiv\sum_{i}\bm{m}_{i}^{\dagger}\bm{m}_{i}$, with $\bm{m}_{i}$ estimated as follows:

[TABLE]

where $\tilde{\Pi}=\bm{\mathcal{Y}}\bm{\mathcal{Y}}^{-1}$ is the projector onto the spherical harmonic space.

The algorithm for the noise covariance estimator proceeds according to the following iterative scheme: Compute $\tilde{\bm{m}_{i}}$ using equation (58), and subsequently $\bm{\mathcal{E}^{-1}}$ to obtain $\widetilde{\mathbf{D}}$ using equation (57), followed by a power spectrum update to obtain $\widetilde{\mathbf{C}}$ via $\widetilde{\mathbf{C}}=\sum_{i}\langle\hat{\bm{m}}_{i}\hat{\bm{m}}_{i}^{\dagger}\rangle/\mathrm{N}_{\mathrm{MC}}$. In the above, we have defined the harmonic representation of the map with $\hat{\bm{m}}_{i}=\bm{\mathcal{Y}}^{-1}\bm{m}_{i}$. We solve equation (57) by implementing fixed point iterations, but this fixed point is not an attractor. We consequently employ the Babylonian method (e.g. Fowler & Robson, 1998; Friberg, 2007) to stabilize the fixed point and obtain an updated $\widetilde{\mathbf{D}}$, as follows:

[TABLE]

We may verify that the fixed point of the above equation is exactly the desired $\mathbf{D}$ matrix. Note that due to the degeneracy between the amplitudes of $\mathbf{C}$ and $\mathbf{D}$, we need to anchor the amplitude of the updated $\widetilde{\mathbf{C}}$ via the required re-scaling.

We therefore solve the above equation (57) iteratively, using $\mathbf{D}^{2}=\sum_{i}(\bm{n}_{i}^{\dagger}\bm{n}_{i})/\mathrm{N}_{\mathrm{MC}}$ as an initial guess. We generate $10^{4}$ MC noise simulations using as template the modulated noise covariance provided by the Planck data analysis pipeline,444Available from http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP_ID=HFI_SkyMap_100_2048_R2.02_full.fits as an estimate for the $\bm{n}_{i}$ above. The map estimates, after only five iterations, for the diagonal and off-diagonal components of the noise covariance matrix, along with their corresponding reference and residual maps, are displayed in Figs. 2 and 3, respectively. Visually, the distinct components of the covariance matrix are adequately recovered, with residuals at the level of $\sim 0.3\%$ and $\sim 6\%$ for the diagonal and off-diagonal components, respectively. As a quantitative diagnostic, we verify the relative deviation in the angular power spectra reconstructed from the maps, as a function of scale, with respect to their reference components: $\sqrt{C_{\ell}(\hat{\mathbf{D}}-\mathbf{D}_{\mathrm{ref}})/C_{\ell}(\mathbf{D}_{\mathrm{ref}})}$, illustrated in Fig. 4. This demonstrates the accuracy of reconstruction of our noise covariance estimator across the range of scales considered, with only five iterations.

5.3 Polarization analysis

In this section, we showcase the application of dante in Wiener filtering polarized CMB maps contaminated with anisotropic correlated noise, and also illustrate its efficacy in generating pure $\mathcal{E}$ and $\mathcal{B}$ maps, guaranteed to be free from any cross-contamination. This corresponds to three distinct runs using the same realization of mock data, generated as described above in Section 5.1, labelled as “WF”, “pure $\mathcal{E}$” and “pure $\mathcal{B}$”, respectively. The “pure $\mathcal{B}$” run yields a “pure” temperature map as by-product, which corresponds to the map of temperature anisotropies without any contribution from the $\mathcal{E}$ modes. We anchor the choice of hyperparameter values, described below, for all three cases.

We make an initial truncation in the power spectrum at $\ell=50$, corresponding to a given value of $\bm{\xi}$ and the algorithm loops through the iterations until the fractional difference between successive iterations has reached a sufficiently low value, at which point $\bm{\xi}$ is reduced by a constant factor according to a given cooling scheme: $\bm{\xi}\rightarrow\bm{\xi}\eta$, where $\eta=2/3$, until $\bm{\xi}\rightarrow\phi\omega\mathbb{1}$. We implement a “weak” criterion for convergence, $\left\|\bm{s}_{i+1}-\bm{s}_{i}\right\|/\left\|\bm{s}_{i}\right\|<\epsilon$, where $\epsilon=10^{-5}$, as a cheap proxy for the strong criterion that is verified a posteriori in Fig. 7.

The reconstructed angular power spectra for the temperature and polarization components are provided in Fig. 5, with the WF solution showing suppressed power on the small scales resulting from the noise and masked regions of the sky, which is a characteristic feature of Wiener filtering. For the temperature anisotropies, depicted in the left panel, the pure $\mathcal{B}$ run yields the temperature power spectra that has been purified with respect to the $\mathcal{E}$ modes, and as such, corresponds to the prediction solely from the temperature data, with no contribution from the polarization component. This pure temperature power spectrum does not display any significant difference compared to the Wiener-filtered one, as expected, due to the relatively small $\mathcal{E}$-mode contribution. The corresponding reconstructed power spectra for the $\mathcal{E}$ modes are depicted in the middle panel of Fig. 5. The pure $\mathcal{E}$-mode power spectrum displays a smooth functional behaviour that matches the shape of the input power spectrum, although substantially suppressed because of the noisy and masked regions and since the discarded temperature contribution is significant. The right panel displays the corresponding reconstructions for the $\mathcal{B}$-mode power spectrum. The pure reconstruction, as in the WF case, suppresses the modes in the low signal-to-noise regime, while also discarding the ambiguous modes, and shows a notable difference on the largest scales.

The real-space maps of the Wiener-filtered, pure $\mathcal{E}$ and $\mathcal{B}$ modes, together with their corresponding simulated maps, are illustrated in the middle and bottom rows of Fig. 6, respectively. The corresponding temperature maps are also displayed in the top row, for completeness. Note that the simulated maps are the generated signal realizations which were subsequently contaminated with anisotropic correlated noise and masked according to Fig. 1. The Wiener-filtered map, as the maximum a posteriori reconstruction, represents the CMB signal content of the data, with the reconstruction of the large-scale modes in the masked areas, based on the information content of the observed sky regions, being a natural consequence of Wiener filtering. Both the Wiener-filtered and pure temperature maps are similar in appearance, as expected from their reconstructed power spectra (cf. left panel of Fig. 5). This is not the case, however, for the reconstructed $\mathcal{E}$ maps, purified with respect to the temperature anisotropies and ambiguous modes. The pure $\mathcal{E}$ map, as a result, has a lower signal amplitude. Both the Wiener-filtered and pure $\mathcal{B}$ maps show the clear reconstruction of the large-scale modes with extremely low amplitude. The striking contrast between the two maps is on the largest scales near the mask, where the pure $\mathcal{B}$ map has lower power. This is consistent with previous work (e.g. Bunn & Wandelt, 2017), where it was found that ambiguous modes have support primarily near the masked regions.

We illustrate the convergence behaviour of the three solutions via the corresponding variations of their residual error given by $\left\|\mathcal{A}\bm{x}-\bm{y}\right\|/\left\|\bm{y}\right\|$, for a linear system of equations given by $\mathcal{A}\bm{x}=\bm{y}$, in Fig. 7. This residual error adequately characterizes the accuracy of the final solution, with the relevant equations as follows:

[TABLE]

and

[TABLE]

with $\bm{x}=\mathcal{S}^{-1/2}\mathfrak{R}\bm{s}$, following the notation from Section 3.2. The coupling matrix $(\bm{\mathcal{Y}}\mathbf{C}\bm{\mathcal{Y}}^{\dagger})^{-1}$ requires Jacobi iterations for accurate evaluation of the above residual error. For the pure $\mathcal{E}/\mathcal{B}$ runs, this is non-trivial as the signal covariance $\mathbf{S}$ should have infinite values, as mentioned at end of Section 4.3. To simplify the residual error evaluations in these cases, we simply set the relevant components of $\mathbf{S}$ to a numerically large value.

As we demonstrated in our previous work, a characteristic feature of the dual messenger algorithm is the monotonic decrease in the residual error as the iterations proceed, as illustrated in the left panel of Fig. 7, thereby demonstrating the unconditional stability of our method. This residual error, as a function of angular scale, for the final solutions from the three runs, are depicted in the middle panel. The Jacobi relaxation schemes employed is required to reduce this error to extremely low values across the range of scales considered.

The $\chi^{2}$ is computed as follows:

[TABLE]

where, as for the residual error evaluations above, we employ Jacobi relaxation for the composite operation $(\bm{\mathcal{Y}}\mathbf{C}\bm{\mathcal{Y}}^{\dagger})^{-1}$. The corresponding $\chi^{2}$ variation for the three different solutions is displayed in the right panel of Fig. 7. In all three cases, the respective $\chi^{2}$ of the dual messenger solutions drop to a final value which matches $\langle\chi^{2}_{\mathrm{d.o.f}}\rangle$, the expectation value of the $\chi^{2}$, given by the number of degrees of freedom (d.o.f), for the final solution. In the absence of masks, $\langle\chi^{2}_{\mathrm{d.o.f}}\rangle$ is given by the total number of harmonic modes of the temperature, $\mathcal{E}$ and $\mathcal{B}$ components. The computation of $\langle\chi^{2}_{\mathrm{d.o.f}}\rangle$ is, however, non-trivial when masks are involved. We estimated $\langle\chi^{2}_{\mathrm{d.o.f}}\rangle$ via Monte Carlo simulations. The convergence diagnostics discussed above, therefore, quantitatively demonstrate the efficacy of dante in performing the three distinct tasks.

Concerning the execution times for the WF, pure $\mathcal{E}$ and pure $\mathcal{B}$ runs, for the specific test case investigated, the algorithm runs to completion on four cores of a standard workstation, Intel Core i5-4690 CPU (3.50 GHz), in around three hours. Note that a conjugate gradient method can, in principle, deal with such anisotropic noise models, provided that an adequate preconditioner can be found and this is the major stumbling block. Devising an appropriate preconditioner for such a complex problem is an extremely challenging task. For instance, the multi-grid preconditioner developed by Smith et al. (2007) at WMAP resolution and sensitivity is already highly non-trivial. The preconditioner-free approach of the dual messenger algorithm, therefore, is the key advantage.

6 Conclusions and outlook

We present a numerically robust and fast code, dante, for pure $\mathcal{E}/\mathcal{B}$ decomposition of CMB polarization maps. It accounts for complex and realistic noise models such as anisotropic correlated noise, encountered in typical CMB experiments such as Planck. dante is an augmented version of our dual messenger algorithm, adapted for the reconstruction of pure full-sky $\mathcal{E}$ and $\mathcal{B}$ maps on the sphere. The algorithm encodes a new method for the pure-ambiguous decomposition, based on a Wiener filtering approach, recently proposed by Bunn & Wandelt (2017), that guarantees no cross-contamination between the two maps. We also developed a noise covariance estimator to reconstruct the components of anisotropic noise covariance from Monte Carlo simulations, as required by the dual messenger algorithm.

We have demonstrated the capabilities of dante in dealing with large data sets and the associated high-dimensional covariance matrices. Moreover, as a preconditioner-free method, it is not hindered by ill-conditioned systems of equations inherent in CMB polarization problems, unlike standard PCG solvers, as demonstrated in KLW18. dante also has an in-built option for drawing constrained Gaussian realizations of the CMB sky, for applications requiring homogeneous coverage of the field of observations. We have not illustrated this particular feature in this work as this was shown previously in KLW18. dante will be rendered public in the near future.

The pure $\mathcal{E}/\mathcal{B}$ decomposition framework implemented in this work, as a maximum a posteriori probability approach, has several advantages over traditional methods. It exploits the sparsity of the $\mathcal{E}/\mathcal{B}$ decomposition in the spherical harmonic basis, rendering the implementation extremely efficient. It is therefore much faster and straightforward than methods relying on the construction of orthonormal bases or wavelet methods that require a certain degree of fine-tuning. Moreover, $\mathcal{E}/\mathcal{B}$ purification in the context of pseudo-$C_{\ell}$ estimators is only feasible when the mask is differentiable up to at least its second derivatives, which is usually achieved via an appropriate apodization (Alonso et al., 2019). An interesting aspect of our approach is that it is not hindered by such limitations.

We have showcased the performance of dante on a realistic mock data set, emulating the features of polarized Planck CMB maps. The next step in this series of investigations is to further augment dante with an adaptive upgrade. Despite the improvements made to render the analysis of high-resolution CMB polarization data sets numerically feasible, the statistically optimal approach for the separation of $\mathcal{E}$ and $\mathcal{B}$ modes requires exact global analyses such as Gibbs sampling. This would, however, require several applications of the Wiener filter to obtain one signal realization conditional on the polarization data (e.g. Larson et al., 2007). The algorithm would therefore benefit from a further level of sophistication. A particularly interesting upgrade is to exploit the hierarchical framework of the dual messenger algorithm by adapting the working resolution progressively during execution, thereby substantially reducing the computation time. We also intend to explore the possibility of employing the dual messenger as a preconditioner in a standard PCG approach, in an attempt to drastically improve the convergence rate. Our algorithm could also be used to provide better examples of Wiener-filtered maps for a filtering based on machine learning, as in Münchmeyer & Smith (2019) with the $J_{1}$ loss function. These examples would be much more expensive than the purely simulation-based approach of the $J_{2}$ loss function, but would also provide a solid validation step.

Ultimately, the underlying objective is to employ this efficient tool in exact global Bayesian analyses of high-resolution and high-sensitivity CMB observations from the latest release of Planck to yield scientific products of significant value and interest. The resulting reconstructed maps may potentially be employed for various applications such as power spectrum reconstruction, estimation of lensing potential, tests for foreground contamination and searches for non-Gaussianity and statistical anisotropy such as hemispherical power anisotropy. Another key aspect is that the features of the real-space pure $\mathcal{B}$ maps allow the characterization of lensing-induced $\mathcal{B}$ modes which go beyond the power spectrum.

dante can be easily applied to other CMB data sets in straightforward fashion without major modifications of the source code, making it a potentially powerful and robust tool for other current and future high-resolution CMB missions such as South Pole Telescope, Advanced ACTPol, Simons Observatory and CMB-S4. The flexibility of the code can nevertheless be exploited in other cosmological contexts, due to the ubiquitous use of the Wiener filter, and even in more general scenarios involving spin field reconstruction.

Acknowledgements

We convey our appreciation to the anonymous reviewer for their remarks and suggestions which helped to improve the overall quality of the manuscript. We thank Franz Elsner and Emory F. Bunn for stimulating discussions. This work has been done within the activities of the Domaine d’Intérêt Majeur (DIM) Astrophysique et Conditions d’Apparition de la Vie (ACAV), and received financial support from Région Ile-de-France. We acknowledge financial support from the ILP LABEX, under reference ANR-10-LABX-63, which is financed by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. GL also acknowledges financial support from the ANR BIG4, under reference ANR-16-CE23-0002. BDW is supported by the Simons Foundation. This work is done within the Aquila Consortium.555https://aquila-consortium.org

Appendix A Spherical harmonic transforms

We provide a brief description of the transformation between pixel and spherical harmonic domain in order to be precise about the notation employed in this work.

Assuming the primary CMB fluctuations to be an isotropic Gaussian random field, the CMB signal can be described as a vector of spherical harmonic coefficients, with the associated signal covariance $\mathbf{S}$ given by:

[TABLE]

where $C_{\ell}$ is the CMB power spectrum. The proper basis to represent isotropic Gaussian random fields on the sphere is described by spherical harmonics. Given a grid on the sphere, i.e. a set of pixel positions $\hat{n}_{p}$, we can transform a field expressed in spherical harmonic (SH) basis, with coefficients $s_{\ell m}$, to one sampled on the sphere via SH synthesis, as follows:

[TABLE]

Formally, the SH synthesis may be expressed as a matrix product, $\bm{s}_{(p)}=\bm{\mathcal{Y}}\bm{a}_{(\ell m)}$, where $\bm{\mathcal{Y}}$ is the synthesis operator that encodes the value of the SHs evaluated at each $\hat{n}_{p}$ of the grid.

Conversely, the transformation from pixel to harmonic basis is referred to as SH analysis, with the analysis operator being an integral related to the synthesis operator:

[TABLE]

It is important to emphasize the scaling operation above and note that the last equalities are valid only for an equal-area pixelization such as HEALPix (Górski et al., 2005). These spherical harmonic transforms are the spherical analogue of Fourier transforms.

Appendix B Modulated correlated noise covariance

We provide a more in-depth derivation of the two dual messenger equations (25) and (26) required for the treatment of modulated correlated noise covariance. The third equation (24) can be derived from equation (23) in straightforward fashion via linear algebraic simplifications.

Equation (5) can be written in its explicit form as:

[TABLE]

with the covariance of the messenger field $\bm{t}$ being $\mathbf{T}=\mathbf{D}(\bm{\mathcal{Y}}\phi\bm{\mathcal{Y}}^{\dagger})\mathbf{D}$. This equation bears a striking resemblance to the standard Wiener filter equation (3) and can therefore be solved via the introduction of an extra messenger field $\bm{v}$ with covariance $\mathbf{V}=\omega(\bm{\mathcal{Y}}\phi\bm{\mathcal{Y}}^{\dagger})\mathbb{1}$, where $\omega\equiv\mathrm{min}(\mathrm{diag}(\mathbf{D}^{2}))$, resulting in the following $\chi^{2}$:

[TABLE]

Minimizing the above $\chi^{2}$ with respect to $\bm{v}$ and $\bm{u}$ yields the following set of equations:

[TABLE]

where, as before, we employ the definition of the coupling matrix $\bm{\mathcal{M}}\equiv\bm{\mathcal{Y}}\bm{\mathcal{Y}}^{\dagger}$. It is more convenient to work with $\tilde{\bm{t}}\equiv\mathbf{D}^{-1}\bm{t}$, such that we can rewrite equation (68) as:

[TABLE]

The preference for the above form is that masked regions do not pose any numerical issue, as $\mathbf{D}^{-1}|_{\mathrm{mask}}=\bm{0}$, such that $\bm{v}|_{\mathrm{mask}}\rightarrow\bm{\mathcal{Y}}\mathbf{B}\bm{u}$.