Wiener filter reloaded: fast signal reconstruction without preconditioning

TL;DR

This paper introduces a dual messenger algorithm for Wiener filtering that efficiently reconstructs signals from large, complex datasets without preconditioning, outperforming existing methods in speed while maintaining accuracy.

Contribution

The paper presents a novel dual messenger algorithm that improves computational efficiency in Wiener filtering for complex data, without requiring preconditioning.

Findings

Matches the performance of the PCG method in reconstruction accuracy.

Runs 2 to 4 times faster than standard messenger and PCG methods.

Successfully applied to simulated CMB temperature data.

Abstract

We present a high performance solution to the Wiener filtering problem via a formulation that is dual to the recently developed messenger technique. This new dual messenger algorithm, like its predecessor, efficiently calculates the Wiener filter solution of large and complex data sets without preconditioning and can account for inhomogeneous noise distributions and arbitrary mask geometries. We demonstrate the capabilities of this scheme in signal reconstruction by applying it on a simulated cosmic microwave background (CMB) temperature data set. The performance of this new method is compared to that of the standard messenger algorithm and the preconditioned conjugate gradient (PCG) approach, using a series of well-known convergence diagnostics and their processing times, for the particular problem under consideration. This variant of the messenger algorithm matches the performance of the PCG method in terms of the effectiveness of reconstruction of the input angular power spectrum and converges smoothly to the final solution. The dual messenger algorithm outperforms the standard messenger and PCG methods in terms of execution time, as it runs to completion around 2 and 3-4 times faster than the respective methods, for the specific problem considered.

Equations

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

$$({\mathbf{S}}^{-1}+{\mathbf{N}}^{-1}){\mathbf{s}}_{\text{WF}}={\mathbf{N}}^{-1}\mathbf{d},$$

$${\chi }_T^2=(\mathbf{d}-\mathbf{t}{)}^{†}{\bar{\mathbf{N}}}^{-1}(\mathbf{d}-\mathbf{t})+(\mathbf{t}-\mathbf{s}{)}^{†}{\mathbf{T}}^{-1}(\mathbf{t}-\mathbf{s})+{\mathbf{s}}^{†}{\mathbf{S}}^{-1}\mathbf{s},$$

$$\left[{\bar{\mathbf{N}}}^{-1}+(\lambda \mathbf{T}{)}^{-1}\right]\mathbf{t}$$

$$\left[{\mathbf{S}}^{-1}+(\lambda \mathbf{T}{)}^{-1}\right]\mathbf{s}$$

$${ϵ}_{i+1}$$

$$=\left[\mathbf{S}(\mathbf{S}+\lambda \mathbf{T}{)}^{-1}\right]\left\{(\mathbf{N}-\mathbf{T}){\left[\mathbf{N}+(\lambda -1)\mathbf{T}\right]}^{-1}\right\}{ϵ}_i.$$

$${\chi }_U^2=(\mathbf{d}-\mathbf{s}{)}^{†}{\mathbf{N}}^{-1}(\mathbf{d}-\mathbf{s})+(\mathbf{s}-\mathbf{u}{)}^{†}{\mathbf{U}}^{-1}(\mathbf{s}-\mathbf{u})+{\mathbf{u}}^{†}{\bar{\mathbf{S}}}^{-1}\mathbf{u},$$

$$\left({\mathbf{N}}^{-1}+{\mathbf{U}}^{-1}\right)\mathbf{s}$$

$$\left({\mathbf{U}}^{-1}+{\bar{\mathbf{S}}}^{-1}\right)\mathbf{u}$$

$${ϵ}_{i+1}$$

$$=[\mathbf{N}(\mathbf{N}+\mathbf{U}{)}^{-1}][(\mathbf{S}-\mathbf{U}){\mathbf{S}}^{-1}]{ϵ}_i.$$

$${\chi }_{\xi }^2=(\mathbf{d}-\mathbf{t}{)}^{†}{\bar{\mathbf{N}}}^{-1}(\mathbf{d}-\mathbf{t})+(\mathbf{t}-\mathbf{s}{)}^{†}{\mathbf{\xi }}^{-1}(\mathbf{t}-\mathbf{s})+{\mathbf{s}}^{†}{\bar{\mathbf{S}}}^{-1}\mathbf{s},$$

$$\left({\bar{\mathbf{N}}}^{-1}+{\mathbf{\xi }}^{-1}\right)\mathbf{t}$$

$$\left({\mathbf{\xi }}^{-1}+{\bar{\mathbf{S}}}^{-1}\right)\mathbf{s}$$

$${\chi }_T^2=(\mathbf{d}-\mathbf{s}{)}^{†}{\left[\mathbf{N}+(\lambda -1)\alpha 1\right]}^{-1}(\mathbf{d}-\mathbf{s})+{\mathbf{s}}^{†}{\mathbf{S}}^{-1}\mathbf{s}.$$

$${\mathbf{s}}_{\text{WF}}(\lambda )$$

$$=\mathbf{S}{\left[\mathbf{S}+\mathbf{N}+(\lambda -1)\alpha 1\right]}^{-1}\mathbf{d},$$

$${\mathbf{s}}_{\text{WF}}(\lambda )-{\mathbf{s}}_{\text{WF}}({\lambda }^{\prime })$$

$$=\mathbf{S}[\mathbf{S}+\mathbf{N}+(\lambda -1)\alpha 1{]}^{-1}\left\{\mathbf{S}+\mathbf{N}+({\lambda }^{\prime }-1)\alpha 1-[\mathbf{S}+\mathbf{N}+(\lambda -1)\alpha 1]\right\}[\mathbf{S}+\mathbf{N}+({\lambda }^{\prime }-1)\alpha 1{]}^{-1}\mathbf{d}$$

$$=\alpha ({\lambda }^{\prime }-\lambda )\mathbf{S}{\left[\mathbf{S}+\mathbf{N}+(\lambda -1)\alpha 1\right]}^{-1}{\left[\mathbf{S}+\mathbf{N}+({\lambda }^{\prime }-1)\alpha 1\right]}^{-1}\mathbf{d}.$$

$${\Vert {\mathbf{s}}_{\text{WF}}(\lambda )\Vert }_{\mathbf{\ell }}=C_{\mathbf{\ell }}(C_{\mathbf{\ell }}+\lambda \alpha {)}^{-1}\Vert \mathbf{d}{\Vert }_{\mathbf{\ell }}.$$

$$\frac{{\Vert {\mathbf{s}}_{\text{WF}}(\lambda )-{\mathbf{s}}_{\text{WF}}({\lambda }^{\prime })\Vert }_{\mathbf{\ell }}}{{\Vert {\mathbf{s}}_{\text{WF}}(\lambda )\Vert }_{\mathbf{\ell }}}$$

$$\le \frac{\alpha |{\lambda }^{\prime }-\lambda |}{C_{\mathbf{\ell }}+{\lambda }^{\prime }\alpha },$$

$$1-ϵ\le \frac{C_{\mathbf{\ell }}}{C_{\mathbf{\ell }}+\lambda \alpha }\left(\frac{\bar{\mathbf{N}}}{\bar{\mathbf{N}}+\lambda \alpha }\right),$$

$$ϵ\le 1-\frac{C_{\mathbf{\ell }}}{C_{\mathbf{\ell }}+\lambda \alpha }=\frac{\lambda \alpha }{C_{\mathbf{\ell }}+\lambda \alpha }.$$

$$\frac{\alpha |{\lambda }^{\prime }-\lambda |}{C_{\mathbf{\ell }}+{\lambda }^{\prime }\alpha }\text{la}\{\begin{array}{c}\frac{\lambda \alpha }{C_{\mathbf{\ell }}+\lambda \alpha }\\ \\ \frac{{\lambda }^{\prime }\alpha }{C_{\mathbf{\ell }}+{\lambda }^{\prime }\alpha }\end{array}$$

$$\left({\mathbf{N}}^{-1}+{\mathbf{U}}_{\mu }^{-1}\right){\mathbf{s}}_{\text{WF}}(\mu )$$

$$\left({\mathbf{U}}_{\mu }^{-1}+{\bar{\mathbf{S}}}_{\mu }^{-1}\right){\mathbf{u}}_{\mu }$$

$${\mathbf{N}}^{-1}\mathbf{d}$$

$$=\left[({\mathbf{N}}^{-1}+{\mathbf{U}}_{\mu }^{-1})({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu })-{\bar{\mathbf{S}}}_{\mu }{\mathbf{U}}_{\mu }^{-1}\right]({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu }{)}^{-1}{\mathbf{s}}_{\text{WF}}(\mu )$$

$$=\left[{\mathbf{N}}^{-1}({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu })+1\right]({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu }{)}^{-1}{\mathbf{s}}_{\text{WF}}(\mu )$$

$$=\left[{\mathbf{N}}^{-1}+({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu }{)}^{-1}\right]{\mathbf{s}}_{\text{WF}}(\mu ),$$

$${\mathbf{s}}_{\text{WF}}(\mu )=\left[{\mathbf{N}}^{-1}+({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu }{)}^{-1}\right]{\mathbf{N}}^{-1}\mathbf{d}.$$

$${\mathbf{s}}_{\text{WF}}(\tilde{\mu })-{\mathbf{s}}_{\text{WF}}(\mu )$$

$$=\left[({\bar{\mathbf{S}}}_{\tilde{\mu }}+{\mathbf{U}}_{\tilde{\mu }})(\mathbf{N}+{\bar{\mathbf{S}}}_{\tilde{\mu }}+{\mathbf{U}}_{\tilde{\mu }}{)}^{-1}-({\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu })(\mathbf{N}+{\bar{\mathbf{S}}}_{\mu }+{\mathbf{U}}_{\mu }{)}^{-1}\right]\mathbf{d}.$$

$${\bar{\mathbf{S}}}_{\tilde{\mu }}={\bar{\mathbf{S}}}_{\mu }+{\Delta }_{\tilde{\mu },\mu }1,$$

Full text

††pubyear: 2017††pagerange: Wiener filter reloaded: fast signal reconstruction without preconditioning–C

1 Introduction

In an era of precision cosmology, the data analysis of state-of-the-art galaxy redshift surveys and CMB experiments with unprecedented levels of sensitivity and resolution poses complex numerical challenges. Fast and robust methods are therefore required to render the data analysis of these large and complex data sets computationally tractable.

One of the most frequently encountered and ubiquitous problems in astrophysics and cosmology (and many other fields in science) is signal reconstruction from noisy data. Solutions to this problem have therefore been researched extensively for the last two centuries (e.g., Gauss, 1809; Jaynes, 1957; Kalman, 1960; Elsner & Wandelt, 2013). The Wiener filter (Wiener, 1949) has emerged as a standard tool for the analysis of large data sets for the inference of high dimensional signals, such as the large-scale structures and CMB problems. It has therefore been employed for large-scale structure analysis problems such as inferring the three dimensional density field from observations (e.g., Zaroubi et al., 1999; Zaroubi, 2002; Erdoǧdu et al., 2004, 2006; Kitaura & Enßlin, 2008; Kitaura et al., 2009; Jasche et al., 2010; Jasche & Lavaux, 2015). In the analysis of CMB data sets, the Wiener filter has been applied to a range of problems, such as the joint inference of temperature fluctuations and power-spectra (e.g., Wandelt et al., 2004; Eriksen et al., 2004; Jewell et al., 2004; O’Dwyer et al., 2004; Smith et al., 2007; Larson et al., 2007; Elsner & Wandelt, 2013; Bunn & Wandelt, 2016).

If we assume a linear model where the data $\bm{d}$ is a combination of the signal $\bm{s}$ and noise $\bm{n}$, i.e.,

[TABLE]

the Wiener filter is defined as the solution to the following equation:

[TABLE]

where $\mathbf{S}$ and $\mathbf{N}$ are the signal and noise covariance matrices, respectively, and $\bm{s}_{\text{\tiny{{WF}}}}$ is the Wiener filter solution. As is evident from Equation (2), the direct numerical implementation of the Wiener filter requires inversion of dense matrices. This task is rendered intractable by the size of modern data sets from state-of-the-art experiments, and therefore represents a computational bottleneck. While it would be extremely convenient if there existed a common basis set, easily accessible by fast transforms, where both $\mathbf{S}$ and $\mathbf{N}$ are sparse, this is often not the case as for instance, the signal and noise covariances may be sparse in Fourier or pixel space, respectively. Some previous approaches relied on the assumption of a homogeneous and isotropic noise distribution to find approximate solutions to the Wiener filter equation (e.g., Hirata et al., 2004; Komatsu et al., 2005; Mangilli & Verde, 2009). Alternative approaches for the exact solution to the Wiener filter equation involve complex numerical algorithms like Krylov space methods, such as preconditioned conjugate gradient (PCG) techniques (eg., Wandelt et al., 2004; Eriksen et al., 2004; Kitaura & Enßlin, 2008, and references therein) for matrix inversion and solving high dimensional systems of linear equations. However, this procedure requires intricate and costly numerical schemes such as preconditioning the linear system with suitable matrices and subsequently demands significant investment in software development. Moreover, this approach involves a further inherent stumbling block as finding efficient preconditioners is a complicated task in itself (Oh et al., 1999) since matrices are often extremely ill-conditioned in typical CMB problems (Eriksen et al., 2004). The inclusion of polarisation data in the analysis further exacerbates this predicament (Larson et al., 2007).

Some previous schemes based on the PCG method involved the use of a combination of block and diagonal preconditioners on large and small angular scales, respectively (Eriksen et al., 2004), or were based on a recursive algorithm, where the conjugate gradient solution on a coarse grid is adopted as the preconditioner on a finer grid (Smith et al., 2007). Although both of these approaches yielded satisfactory performance with the analysis of WMAP data (e.g., Eriksen et al., 2008), they were found to be too computationally intensive for the high resolution and sensitivity analysis of Planck data. Seljebotn et al. (2014) recently proposed a multi-level solver for Gaussian constrained realisations of the CMB, which is fast once a suitable preconditioner is found but requires careful tuning and costly precomputations in terms of both computing power and memory requirements. This approach is therefore less attractive when we have to solve many different systems, each requiring a specific preconditioner.

Elsner & Wandelt (2013) presented an alternative approach by devising the messenger algorithm which bypasses the need for a preconditioner, as described in the following section. The evaluation of the Wiener filter from Equation (2) via both the standard and a new formulation that is dual to the messenger algorithm constitutes the crux of this work.

The paper is organised as follows. In Section 2, the underlying principles of the standard messenger method are outlined, followed by a description of the new dual messenger algorithm. We then test our new technique on an artificially generated CMB data set in Section 3, and follow up by investigating its performance in terms of convergence, computation time and stability, and draw comparisons to the standard messenger scheme and the popular PCG method in Section 4. Finally, in Section 5, we summarise the main aspects of our findings and discuss the areas of applications where the potential of our new algorithm can be fully exploited. In Appendices A and B, we illustrate the rationale behind the schemes implemented in this work for fast convergence, followed by a brief review of the PCG method in Appendix C.

2 The Messenger algorithms

2.1 The standard messenger algorithm

Elsner & Wandelt (2013) proposed a high precision, iterative algorithm for the solution to the full Wiener filter equation, while being numerically efficient and straightforward to implement. Conceptually, the key idea is to introduce a stochastic auxiliary field $\bm{t}$, the so-called “messenger" field, with covariance $\mathbf{T}$, where $\mathbf{T}$ is proportional to the identity matrix. Taking advantage of the useful property of the identity matrix being invariant under any orthogonal basis transformation, the messenger field acts as an intermediate in transformations between different preferred orthogonal bases, in which signal and noise covariance matrices are expressed conveniently, i.e., are sparse. As a result, although directly applying combinations like $(\mathbf{S}+\mathbf{N})^{-1}$ may not be possible, we can always apply expressions like $(\mathbf{S}+\mathbf{T})^{-1}$ and $(\mathbf{N}+\mathbf{T})^{-1}$, irrespective of the basis chosen to render $\mathbf{S}$ and $\mathbf{N}$ sparse. Under such a scheme, the information from the data is transmitted to the signal via the messenger field that can be transformed efficiently from one basis representation to another, thereby obviating the requirement to apply the inverse Wiener covariance matrix to data.

With the introduction of the messenger field $\bm{t}$, the modified $\chi^{2}$, where the posterior probability distribution of $\bm{s}$ is proportional to $\exp\left({-\chi^{2}/2}\right)$, is as follows:

[TABLE]

where we defined $\bar{\mathbf{N}}\equiv\mathbf{N}-\mathbf{T}$, and we choose the covariance matrix of the auxiliary field $\bm{t}$ according to $\mathbf{T}=\alpha\mathbb{1}$, where $\alpha\equiv\textup{min}(\textup{diag}(\mathbf{N}))$. Minimising with respect to $\bm{s}$ and $\bm{t}$ leads to the following two equations:

[TABLE]

where we also introduced a scalar parameter $\lambda$ whose purpose is to accelerate convergence. In the limit of $\lambda=1$, the above system of Equations (4) and (5) reduces to the usual Wiener filter Equation (2), as shown in Appendix A. The messenger algorithm basically involves the following steps, as outlined in Algorithm 1. We initialise the vectors $\bm{s}$ and $\bm{t}$ with zeros, and choose an initial high value of $\lambda$. We first solve Equation (4) for the messenger field $\bm{t}$ in the basis defined by $\mathbf{N}$, the noise covariance matrix. Then, we change to a basis where the signal covariance matrix $\mathbf{S}$ has a sparse representation, such as Fourier space. Next, we solve for $\bm{s}$ using Equation (5) and the resulting $\bm{t}$ from the previous step. Finally, we transform the result back to the original basis. The signal reconstruction converges to the Wiener filter solution, i.e., $\bm{s}\rightarrow\bm{s}_{\text{\tiny{{WF}}}}$, as $\lambda\rightarrow 1$.

If $\epsilon_{i}$ is the residual at the $i$th step, then at the following $(i+1)$th step, the corresponding residual is

[TABLE]

For the case of homogeneous noise, i.e., $\mathbf{N}\propto\mathbb{1}$, the system is solved exactly in a single step. A key observation is that $\left|\epsilon_{i+1}\right|<\left|\epsilon_{i}\right|$ for all $i$ since the terms in brackets are less than unity, resulting in unconditional convergence of the signal reconstruction $\bm{s}$ to the Wiener filter solution $\bm{s}_{\text{\tiny{{WF}}}}$. A careful inspection of Equation (7) yields the following insight: Convergence is fast for low noise pixels and modes with low signal prior variance while the converse is also true, i.e., the system converges slowly for high prior variance and high noise pixels. In the latter regimes, a high value of $\lambda$ would speed up convergence. Elsner & Wandelt (2013) have shown that it is possible to find a cooling scheme for $\lambda$, where $\lambda\gg 1$ initially, to smoothly bring the algorithm to the final solution ($\lambda=1$). Here, we reduce $\lambda$ by a constant factor $(1/\eta)$, where $1/2\leq\eta\leq 1$. The rationale behind the cooling scheme adopted in this work is laid out in Appendix A.

The standard messenger method has enjoyed considerable success, thereby establishing its credentials as a reliable method for fast Wiener filtering, without having recourse to a preconditioner. Essentially, it speeds up the computation of the Wiener filter when the noise is strongly inhomogeneous and/or the data is masked. Elsner & Wandelt (2012); Elsner & Wandelt (2013) applied the messenger algorithm on CMB data from WMAP satellite and found the final map to be accurate to about 1 part in $10^{5}$, compared to standard conjugate gradient solvers. Even with the inclusion of polarisation data in the analysis, the messenger algorithm maintained its efficiency. Mangilli et al. (2013) implemented the messenger algorithm to produce Wiener-filtered simulations of non-Gaussian CMB maps with the lensing-integrated Sachs Wolf bispectrum signal. The messenger method was also employed in Planck data analysis, specifically for inverse covariance filtering of CMB maps at high angular resolutions (Planck Collaboration et al., 2014a, b). Jasche & Lavaux (2015) implemented a Gibbs sampling adaptation of the messenger algorithm in a simple, easy to implement but efficient algorithm for Bayesian large-scale structure inference, specifically aiming at the joint inference of cosmological density fields and power spectra for linear data models. Anderes et al. (2015) adopted a similar approach in their Bayesian hierarchical modelling of the CMB gravitational lensing. Alsing et al. (2016a) also implemented the Gibbs-messenger sampling adaptation developed by Jasche & Lavaux (2015) for Bayesian hierarchical modelling of cosmic shear power spectrum inference and eventually for cosmological parameter inference (Alsing et al., 2016b). The works of Jasche & Lavaux (2015), Anderes et al. (2015) and Alsing et al. (2016a) show that the messenger algorithm can be successfully adapted for high resolution conditional Gaussian sampling using Markov Chain Monte Carlo techniques, demonstrating the flexibility of the method.

2.2 The dual messenger algorithm

In the above messenger algorithm, we made use of a relation to trivially split the noise into two components: one with a trivial covariance matrix and the other as a fluctuating component over the sky. An alternative approach is to introduce the auxiliary field at the level of the signal in a complementary formalism to the standard messenger framework. Due to the two schemes being complementary to each other, the new algorithm is referred to as the dual messenger algorithm. The corresponding log-posterior to be optimised then becomes

[TABLE]

where, analogous to the standard approach, $\mathbf{U}=\nu\mathbb{1}$ with $\nu\equiv\mathrm{min}(\textup{diag}(\mathbf{S}))$, and the covariance of the auxiliary field, $\bar{\mathbf{S}}\equiv\mathbf{S}-\mathbf{U}$. We derive the corresponding equations that must be satisfied by $\bm{s}$ and $\bm{u}$ at the minimum of ${\chi^{2}_{\tiny{U}}}$:

[TABLE]

The dual messenger algorithm has interesting convergence properties. The amount of reduction of the residual at each iteration is given by

[TABLE]

This provides the basis for the following mechanism: We artificially truncate the spectrum $\mathbf{S}$ to some lower initial value of $\ell_{\mathrm{iter}}$ that corresponds to a covariance $\mu$ and bring $\ell_{\mathrm{iter}}$ slowly to $\ell_{\mathrm{max}}$ corresponding to our final covariance $\nu$. So, essentially, we vary the covariance $\mathbf{U}$ via a cooling scheme to bring $\mu\rightarrow\nu$, where, in the limit $\mu=\nu$, the above system of Equations (9) and (10) reduces to the usual Wiener filter Equation (2), as shown in Appendix B. This results in a redefinition of $\bar{\mathbf{S}}$ using the Heaviside function as $\bar{\mathbf{S}}=\Theta(\mathbf{S}-\mathbf{U})$, as described quantitatively in Appendix B (cf. Equations (32) and (33)). For $\ell\la\ell_{\mathrm{iter}}$, the ratio given by Equation (12) is always convergent. The algorithm consists of the following steps. We initialise the vectors $\bm{s}$ and $\bm{u}$ with zeros. We begin iterations with an initial value of $\ell_{\mathrm{iter}}$ and correspondingly $\mu$, and iterate until $\mu\rightarrow\nu$ at $\ell_{\mathrm{iter}}=\ell_{\mathrm{max}}$, i.e., $\bm{s}\rightarrow\bm{s}_{\text{\tiny{{WF}}}}$, in accordance with a chosen convergence criterion for each value of $\mu$. In analogy with the messenger technique, we need to adopt a cooling scheme for $\mu$. The idea is to reduce $\mu$ by a constant factor $(1/\beta)$, where $0<\beta<1$, which matches the convergence speed of one iteration with the modification in $\mu$. The rationale behind the cooling scheme tailored for the dual messenger algorithm is illustrated in Appendix B.

However, numerically, the above algorithm does not result in the correct final solution due to the continuous mode of the signal, i.e., the zero eigenvalue in the signal covariance $\mathbf{S}$. We therefore require $\mu\rightarrow\nu=0$ to obtain the proper solution at the end, which cannot be accommodated by the dual messenger scheme above. To remedy this numerical predicament, we introduce an extra degree of freedom, $\alpha$, in the system, where $\alpha\equiv\textup{min}(\textup{diag}(\mathbf{N}))$, thereby incorporating aspects of the standard messenger method into the dual messenger scheme. This leads to the following $\chi^{2}$:

[TABLE]

where $\bm{\xi}=(\alpha+\mu)\mathbb{1}=\xi\mathbb{1}=\mathbf{T}+\mathbf{U}$, with $\mu$, $\mathbf{T}$, $\bar{\mathbf{S}}$ and $\bar{\mathbf{N}}$ inheriting their previous definitions. The corresponding set of equations to be solved iteratively is then:

[TABLE]

If $\alpha=0$, we recover the usual dual messenger scheme, while setting $\mu=0$ yields the standard messenger algorithm. The definition of $\bm{\xi}$ implies that the cooling scheme described above still applies to this hybrid method. As outlined in Algorithm 2, we proceed in similar steps as described above for the previous scheme, except that here we reduce the norm of $\bm{\xi}$ by the factor ($1/\beta$) and iterate until $\xi\rightarrow\alpha$, at which point $\mu=0$ and we obtain the proper solution as desired.

3 Application to Cosmic Microwave background

We generated artificial CMB data, in a realistic scenario, by drawing Gaussian random fields on a 2d flat sky with grid resolution, $N_{\textup{pixels}}=512^{2}$, and angular extent, $L=10.0$ degrees. This true simulated map is contaminated by Gaussian white noise with covariance matrix, $\mathbf{N}=(64.0\>\mu\mathrm{K^{2}})\>\mathbb{1}$, via a linear model as given in Equation (1). We made use of CAMB111http://camb.info (Lewis et al., 2000) to generate the input angular power spectrum from which the CMB signals are drawn. We assume a standard $\Lambda$CDM cosmology with the set of cosmological parameters ($\Omega_{m}=0.32$, $\Omega_{\Lambda}=0.69$, $\Omega_{b}=0.05$, $h=0.67$, $\sigma_{8}=0.83$, $n_{s}=0.97$) from Planck (Planck Collaboration et al., 2016). A portion of the flat sky is masked by imposing a central square patch of angular extent $5.0$ degrees with a noise covariance of $10^{6}\mathbf{N}$. The three different methods, messenger, dual messenger and PCG, are applied to this set-up to investigate their effectiveness and efficiency of signal reconstruction. We implement the same “weak” Cauchy convergence criterion, $\left\|\bm{s}_{i+1}-\bm{s}_{i}\right\|/\left\|\bm{s}_{i}\right\|<\epsilon$, where $\epsilon=10^{-6}$, in all the algorithms, to ensure unbiased results. As a reference, we make use of the PCG method with the more stringent $\epsilon=10^{-9}$ to provide results against which the other methods can be compared. Currently regarded as the standard Wiener filtering technique by the scientific community, the PCG method is the natural choice for providing a reference solution in the absence of an exact analytic solution. By imposing a more stringent convergence criterion, we ensure that the PCG solution is close to the ground truth.

The true CMB simulated signal and the reconstructed map obtained using the dual messenger algorithm are displayed in Figure 1. The characteristic feature of Wiener filtering, i.e., the extrapolation of signal into the mask, can be distinctly observed from the reconstructed map. At an initial cursory glance, the reconstructed maps from all three methods appear rather similar, and hence, only the dual messenger reconstruction is displayed in Figure 1. This observation holds for the low noise regions but under scrutiny, the reconstruction in masked regions shows some slight differences. The residual maps depicted in Figure 3, generated by computing the difference, over the full sky, between our reference map and the corresponding reconstructed map obtained using each method, help to discern these differences. The messenger technique provides the most accurate reconstruction in the masked region as illustrated by its residual map, with less than $3\%$ residuals. The dual messenger and PCG reconstructions yield around $4\%$ and $9\%$ residuals, respectively, with the residuals obtained via the former algorithm being concentrated on the edges of the mask.

The effectiveness of reconstruction of the dual messenger technique on both small and large scales is also manifest from Figure 2 which shows the reconstructed power spectra from the different algorithms. The corresponding power spectra recovered via the messenger, dual messenger and PCG methods are all in good agreement, on all scales, with the reference one, computed using the PCG scheme with the more stringent convergence criterion of $\epsilon=10^{-9}$. The characteristic feature of Wiener filtering where the recovered power on the small scales in the low signal to noise regime is suppressed is also observed from Figure 2. Standard Wiener filtering algorithms usually encounter difficulties in dealing with masked regions having infinite noise. By masking $25\%$ of the simulated map with noise of the order $10^{6}$ higher than the unmasked region, we are investigating the worst case performance aspect of the dual messenger algorithm. Hence, we do not expect the algorithm to encounter any difficulties in Wiener filtering maps with high levels of inhomogeneous noise.

The dual messenger algorithm is therefore relevant for current and future high resolution CMB experiments such as South Pole Telescope, Advanced ACTPol, Simons Observatory and CMB-S4. The application of the algorithm can be extended to problems involving the polarisation of the CMB. The formalism remains unchanged for spin field reconstruction, although the numerical implementation is less trivial due to the correlation between the temperature and polarisation components of the signal covariance. The Wiener filtering of polarised CMB data with more complex noise models will be subjected to future investigation to further showcase the efficiency of the dual messenger algorithm in treating complex CMB problems.

4 Numerical Convergence and Stability

We now investigate the convergence properties of the dual messenger algorithm and also provide a more in depth study of the performance of the standard messenger scheme. We quote the usual statistics for convergence and the change in $\chi^{2}$ of the posterior probability density between successive iterations for each method, so that unbiased comparisons of their efficiency and effectiveness can be drawn.

Figure 4 shows the residual error at each iteration for the different methods, as a function of the number of iterations. For the messenger algorithm, we relax the convergence criterion $\epsilon$ for high values of $\lambda$ and reduce $\epsilon$ to $10^{-6}$ as $\lambda\rightarrow 1$ in a three-step procedure, as depicted in Figure 4 in dashed red lines. Here, we adopt a cooling scheme where $\lambda$ is reduced by a constant factor of $4/3$, i.e, $\eta=3/4$, until $\lambda=1$. This consequently ensures that the $\chi^{2}$ decreases rapidly, thereby bringing us closer to the final solution with a smaller number of iterations, resulting in faster convergence.

For the dual messenger algorithm, we implemented a cooling scheme for $\xi$ with $\beta=3/4$, which results in fast convergence while providing accurate results. So, we pick an initial value of $\mu$ corresponding to an initial value of $\ell_{\mathrm{iter}}$, and therefore $\xi$, and iterate until convergence, as dictated by a given Cauchy criterion. We repeat this iterative procedure, in accordance with the aforementioned cooling scheme, until $\xi=\alpha$, at which point we have the Wiener filter solution. Again, we impose less stringent convergence criterion at higher values of $\xi$, as shown in Figure 4 in dashed blue lines.

Figure 4 essentially illustrates the convergence of the different methods. The dual messenger algorithm requires slightly fewer iterations to converge to the final solution than the PCG approach, while the messenger technique requires nearly twice as many iterations as its dual counterpart for convergence. In terms of wall-clock times, the dual messenger has the definite upper hand, as it runs to completion in around 258 seconds whereas the messenger and PCG methods have corresponding wall-clock times of roughly 435 seconds and 663 seconds, i.e., reconstruction via the dual messenger algorithm is nearly three and two times faster than the PCG and messenger methods, respectively. All computations were run on a single core of an Intel Core i5-4690 CPU (3.50 GHz). Both messenger algorithms possess the same algorithmic complexity and memory requirements. They require two Fourier transforms, $\mathcal{O}(N_{\textup{pixels}}\log N_{\textup{pixels}})$, and two scalar multiplications corresponding to algebraic operations of $\mathcal{O}(N_{\textup{pixels}})$, per iteration. In terms of memory requirements, two vectors of size $N_{\textup{pixels}}$ must be temporarily stored in memory. In comparison, the PCG method requires three Fourier transforms and ten scalar multiplications per iteration, and temporary storage of eight vectors of dimension $N_{\textup{pixels}}$ in memory.

We also performed a consistency check by verifying the variation of the residual error given by $\left\|\mathbf{\mathcal{A}}\bm{x}-\bm{y}\right\|/\left\|\bm{y}\right\|$. This is usually adopted as a convergence criterion for PCG computations. The variation of this residual error, analogous to Figure 4, is illustrated in Figure 5. The corresponding residual errors obtained via the different algorithms all drop below their respective Cauchy convergence thresholds implemented. The oscillatory behaviour of the PCG solution, displayed in Figures 4 and 5, is mainly due to the re-initialisation step after every two hundred iterations in the algorithm, as described in Appendix C. However, there are also some oscillations in the residual errors due to the PCG method being susceptible to instabilities sourced by numerical noise. The oscillations in the residual errors in Figure 4 of the two messenger solutions are however due to their respective cooling schemes, resulting in transitions in the systems of equations with the varying covariances of the auxiliary fields (cf. Equations (20) and (37)), with the peaks produced coinciding with these transitions. It is important to note that the residual errors always drop sharply after the peaks, thereby demonstrating the unconditional stability of the messenger algorithms. From numerical experiments, the messenger techniques have proven to be far more stable than the PCG method for nearly degenerate systems.

Due to the Wiener filter being the maximum a posteriori solution, the $\chi^{2}$ of the intermediate solution can be regarded as a useful convergence diagnostic. The $\chi^{2}$ variation as a function of number of iterations is displayed in Figure 6, with the left and right panels correspondingly showing the convergence of the messenger and dual messenger algorithms. The cooling scheme for $\lambda$ implemented in the standard messenger technique causes the $\chi^{2}_{\mathrm{M}}$ to drop rapidly with each change in $\lambda$. Intuitively, this decrease of $\lambda$ via a series of such steps seems reasonable since the $\chi^{2}$ decreases sharply with each change in $\lambda$, and then reaches a plateau for a given $\lambda$ until the latter decreases further, as evidenced in Figure 2 of Elsner & Wandelt (2012). The $\chi^{2}_{\mathrm{DM}}$ of the dual messenger algorithm has a similar behaviour. The $\chi^{2}$ in both cases finally attains the $\chi^{2}_{\mathrm{ref}}$ of the reference PCG method with $\epsilon=10^{-9}$, with $\Delta\chi^{2}_{\mathrm{M}}/\chi^{2}_{\mathrm{ref}}=9.2\times 10^{-6}$ and $\Delta\chi^{2}_{\mathrm{DM}}/\chi^{2}_{\mathrm{ref}}=9.0\times 10^{-6}$, where $\chi^{2}_{\mathrm{ref}}=4.2\times 10^{4}$.

Another important convergence diagnostic is the variation of the relative error, $C_{\ell}(\bm{s}_{\text{\tiny{{WF}}}}-\bm{s}_{\text{\tiny{{WF}}}}^{\textrm{ref}})/C_{\ell}(\bm{s}_{\text{\tiny{{WF}}}}^{\textrm{ref}})$, computed over the full sky, as a function of scale, $\ell$, depicted in Figures 7, 8 and 9 for the messenger, dual messenger and PCG algorithms, respectively. The relative errors on the small scales are of the order $10^{-12}$, $10^{-12}$ and $10^{-15}$, and conversely, on the large scales, $10^{-3}$, $10^{-3}$ and $10^{-2}$, correspondingly, for the messenger, dual messenger and PCG methods. The messenger technique displays smooth and nearly uniform convergence on small and intermediate scales, with the relative error dropping below $10^{-6}$, while remaining below $10^{-3}$ for the largest scales. For the dual messenger scheme, the corresponding convergence rate highlights the hierarchical fashion in which the solution is computed, while yielding similar final relative error across all scales as the messenger algorithm. The PCG method has the lowest relative error on the smallest scales, although this may be biased by the fact that the reference method is also a PCG, but remains inferior to both messenger methods on the largest scales. This is consistent with the significant residuals resulting from the PCG reconstruction, as observed in the previous section (cf. Figure 3). We stress that all relative errors above are computed with respect to the reference PCG method with $\epsilon=10^{-9}$.

We also carry out a series of additional runs with various values of the convergence criterion $\epsilon$ to investigate the wall-clock times required for convergence in each case, thereby providing a more in-depth picture of the performance of the different algorithms. As illustrated in Figure 10, the dual messenger technique converges faster than the other methods, except for the extreme values of $\epsilon$. It is also interesting to note that we can further reduce its execution time by lowering the factor $\beta$ for the cooling scheme without degrading the accuracy of results significantly. For instance, choosing $\beta=1/2$ reduces the number of iterations required, and therefore computation time, by around $25\%$ for the case $\epsilon=10^{-6}$.

5 Conclusions

We presented a new formulation that is dual to the recently developed messenger algorithm, where, unlike in the standard approach, the auxiliary field is introduced at the level of the signal and is consequently associated with the signal reconstruction instead of the noise. This new iterative solver provides another pathway to solve the ill-posed problem of Wiener filtering, frequently encountered in several applications in cosmology and astrophysics.

We tested our new method on a simulated CMB data set and its performance was evaluated in terms of effectiveness of reconstruction, convergence properties, processing time and stability. The dual messenger scheme is shown to match the accuracy of reconstruction of the standard messenger and PCG methods on all scales. Regarding the convergence of the algorithm, it is shown to perform smoothly over all scales, with the relative error being sufficiently low even on the largest scales. For the specific problem under consideration, the dual messenger algorithm has a definite edge over the standard messenger scheme and the PCG approach in terms of computation time. The efficiency and effectiveness of this new technique in calculating the Wiener filter solution of general data sets has therefore been demonstrated. We also provided further insight into the mechanism of the standard messenger method so that it can be optimised for data analysis by the scientific community.

The dual messenger algorithm, like its predecessor, does not require an ingenious choice of preconditioner and is straightforward to implement, robust and flexible. It is also capable of taking into account inhomogeneous noise distributions and arbitrary mask geometries. A key aspect of this technique is that it computes the Wiener filter solution in a hierarchical manner due to the thresholding of the signal covariance matrix inherent in the algorithm. Given the success of the standard messenger method, this new dual messenger technique is highly promising and may therefore be adapted for high precision large-scale data analysis methods in cosmology and astrophysics.

Other performance improvements to this algorithm can be obtained by adapting the working resolution such that the Nyquist frequency is always slightly higher than the current $\ell_{\mathrm{iter}}$ considered in the dual messenger method. This would consequently reduce the number of operations required for Fourier transform. We are also currently working on a generalisation of both the messenger and dual messenger algorithms in a combined approach that may further optimise the performance and efficiency compared to traditional techniques of solving the Wiener filter problem.

Acknowledgements

The authors thank the anonymous reviewer for his constructive comments which helped to improve the manuscript. The authors also thank Franz Elsner for his comments on a draft version of the paper. The authors 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 “Programme National de Cosmologie et Galaxies” (PNCG) of CNRS/INSU, France, and INPHYNITI programme of Interdisciplinary division of CNRS.

Appendix A Cooling Scheme for Messenger Algorithm

In this section, we illustrate the rationale behind the cooling scheme for the scalar parameter $\lambda$ adopted for the standard messenger technique. In the computations below, we assume invertibility of matrices everywhere. Using $\mathbf{T}=\alpha\mathbb{1}$, where $\alpha=\textup{min}(\textup{diag}(\mathbf{N}))$, the ${\chi^{2}_{\tiny{T}}}$ for the messenger scheme from Equation (3) can be written as

[TABLE]

Hence, in the messenger approach, we have $\mathbf{N}\rightarrow[\mathbf{N}+(\lambda-1)\alpha\mathbb{1}]$ and plugging this into the Wiener filter given by Equation (2) yields

[TABLE]

where we assume that all matrices are invertible. As a consistency check, for $\lambda=1$, Equation (17) reduces to Equation (2), the usual Wiener filter equation. The difference between the solutions at two consecutive values of $\lambda$ is given by

[TABLE]

Now, we consider limiting cases of Equation (18). We assume a homogeneous noise distribution, i.e., $\mathbf{N}=\alpha\mathbb{1}$. From Equation (17), we have

[TABLE]

Using Equations (18) and (19), the relative error can then be computed as follows:

[TABLE]

where, for an arbitrary matrix $\mathbf{M}$, $\left\|\mathbf{M}\right\|_{\bm{\ell}}$ is the subset of $\left\|\mathbf{M}\right\|$ over subspace $\bm{\ell}$. From Equation (6), we have the fractional error reduction, after one iteration, given by

[TABLE]

where the term in parentheses will always favour rapid convergence. The other term is strongly convergent on small scales, but to give bounds to the maximum convergence speed of “slow modes”, we focus on the latter, leading to

[TABLE]

To favour convergence, from Equations (20) and (22), we require the two conditions:

[TABLE]

so that the error reduction in one iteration matches the change in solution arising from two consecutive values of $\lambda$. This results in the following two constraints: $|\lambda^{\prime}-\lambda|\leq\lambda^{\prime}$ and $|\lambda^{\prime}-\lambda|\leq\lambda$. The former is automatically satisfied, while the latter leads to the following bounds: $\lambda^{\prime}/2\leq\lambda\leq\lambda^{\prime}$. Therefore, for the cooling scheme, we can write $\lambda=\eta\lambda^{\prime}$, where $1/2\leq\eta\leq 1$, to improve convergence on all scales. This is the motivation behind our choice of $\eta=3/4$ in this work.

Appendix B Cooling Scheme for Dual messenger Algorithm

As in the previous section, we wish to analytically determine the cooling scheme that would improve convergence for the dual messenger algorithm. We first derive the Wiener filter solution for the dual messenger, i.e., the analogue of Equation (17) for this scheme. We begin by writing down Equations (9) and (10), showing the dependence on the power spectrum truncation $\mu$:

[TABLE]

Solving for $\bm{s}_{\text{\tiny{{WF}}}}(\mu)$ via the following steps:

[TABLE]

so that finally, we have

[TABLE]

In the limit $\mu\rightarrow\nu$, where $\nu\equiv\mathrm{min}(\textup{diag}(\mathbf{S}))$, $\bar{\mathbf{S}}+\mathbf{U}\rightarrow\mathbf{S}$, reducing Equation (27) to the usual Wiener filter Equation (2), implying consistency. When we change the truncation of the spectrum, $\mu\rightarrow\tilde{\mu}$, for $\tilde{\mu}<\mu$, we obtain a new solution $\bm{s}_{\text{\tiny{{WF}}}}(\tilde{\mu})$. The counterpart of Equation (20) is then

[TABLE]

Here, we consider Fourier space, so that $\mathbf{U}_{\tilde{\mu}}\rightarrow\tilde{\mu}$ and $\mathbf{U}_{\mu}\rightarrow\mu$. We can write the truncated signal covariance matrix as

[TABLE]

where $\Delta_{\tilde{\mu},\mu}=\tilde{\mu}-\mu$ is the portion of the power spectrum bounded by $\tilde{\mu}$ and $\mu$, while the corresponding truncated signal covariances can be represented by Heaviside functions as follows:

[TABLE]

where, for a matrix $\mathbf{M}=\mathbf{P}\Lambda\mathbf{P}^{-1}$, after applying a basis transformation, with $\Lambda$ being diagonal,

[TABLE]

and

[TABLE]

Using Equation (29) in Equation (28) results in

[TABLE]

Using Equation (29) and substituting the following form of Equation (27),

[TABLE]

in Equation (34) leads to

[TABLE]

such that finally we obtain the following equation for the relative error:

[TABLE]

Now, for $\mu\rightarrow\tilde{\mu}$, there are three distinct regimes of relevance and we consider each case below. To investigate the convergence behaviour in the different regimes, we make use of the following equation, obtained by plugging Equations (30) and (31) in Equation (29),

[TABLE]

The three regimes are:

• •

$\mathbf{S}_{\bm{\ell}}>\mu\mathbb{1}$,

where $\Theta(\mathbf{S}-\tilde{\mu})=\mathbf{S}_{\bm{\ell}}-\tilde{\mu}\mathbb{1}$, and $\Theta(\mathbf{S}-\mu)=\mathbf{S}_{\bm{\ell}}-\mu\mathbb{1}$, so that $\Delta_{\tilde{\mu},\mu}=\mu-\tilde{\mu}$, and this causes the second term in Equation (37) to vanish, implying that the relative error is not affected by our choice of $\tilde{\mu}$. This shows that the Wiener filter solution is naturally computed in a hierarchical fashion using the dual messenger algorithm.

• •

$\tilde{\mu}\mathbb{1}<\mathbf{S}_{\bm{\ell}}<\mu\mathbb{1}$,

where $\Theta(\mathbf{S}-\tilde{\mu})=\mathbf{S}_{\bm{\ell}}-\tilde{\mu}\mathbb{1}$, and $\Theta(\mathbf{S}-\mu)=0$, so that $\Delta_{\tilde{\mu},\mu}=\mathbf{S}_{\bm{\ell}}-\mu\mathbb{1}$, leading to the following relative error:

[TABLE]

The first term behaves as a constant, i.e., $\left\|\mathbf{N}(\mathbf{N}+\bar{\mathbf{S}}_{\mu}+\mu)^{-1}\right\|\sim\alpha^{\prime}$, so that the relative error can be approximated as

[TABLE]

since $\bar{\mathbf{S}}_{\tilde{\mu}}+\tilde{\mu}\sim C_{\bm{\ell}}$. Also, $\alpha^{\prime}\sim 1$, and to favour convergence, we want this change to be as large as possible but sufficiently small such that iterating the solution results in rapid decay of the error, i.e., the change in the solution due to changing $\mu$ should be matched to the change when iterating the solution. If $\mu=\beta C_{\bm{\ell}}$, choosing $0<\beta<1$ would therefore improve convergence by avoiding the early freeze of modes in the iteration scheme. This served as the basis for our choice of $\beta=3/4$ in this work.

• •

$\mathbf{S}_{\bm{\ell}}<\tilde{\mu}\mathbb{1}$,

where $\Theta(\mathbf{S}-\tilde{\mu})=0=\Theta(\mathbf{S}-\mu)$, so that $\Delta_{\tilde{\mu},\mu}=0$, and again using the approximation $\bar{\mathbf{S}}_{\tilde{\mu}}+\tilde{\mu}\sim C_{\bm{\ell}}$ yields

[TABLE]

So, the overall convergence behaviour can be quantitatively described by:

[TABLE]

and this leads to the following interpretation: For higher values of $\mu$ on large scales, the third term in Equation (42) dominates since $\left|\tilde{\mu}-\mu\right|$ is large. But for the final truncations, $\tilde{\mu}\sim\mu$, so this regime is saturated, while on small scales, the second term dominates as $\mu\gg C_{\bm{\ell}}$, so the relative error continues to drop till the final truncation. The cooling scheme described above applies naturally to the hybrid version of the dual messenger method.

Appendix C Brief Review of Preconditioned Conjugate Gradient Method

In general, the PCG approach consists of solving the following set of linear equations:

[TABLE]

where $\mathbf{\mathcal{A}}$ is usually a very large matrix. We wish to avoid computing the inverse of this dense matrix and to this end, in the traditional PCG scheme, we make use of a sparse matrix known as the preconditioner $\mathbf{\mathcal{M}}$, which is the approximate inverse of $\mathbf{\mathcal{A}}$, i.e., $\mathbf{\mathcal{M}}\approx\mathbf{\mathcal{A}}^{-1}$, as follows:

[TABLE]

The Wiener filter Equation (2) can be rewritten as

[TABLE]

such that $\mathbf{\mathcal{A}}=1+\mathbf{S}^{\frac{1}{2}}\mathbf{N}^{-1}\mathbf{S}^{\frac{1}{2}}$, $\bm{x}=\mathbf{S}^{-\frac{1}{2}}\bm{s}_{\text{\tiny{{WF}}}}$ and $\mathbf{y}=\mathbf{S}^{\frac{1}{2}}\mathbf{N}^{-1}\bm{d}$, in accordance with Equation (43). For the case considered in this work, we make use of a suitable preconditioner, which is diagonal in Fourier space, as follows:

[TABLE]

where $n_{\bm{k}}$ are the eigenvalues of the inverse noise covariance matrix, $\mathbf{N}^{-1}$. We then implement this preconditioner in a PCG algorithm (e.g., Golub & Van Loan, 1996). We make an initial guess $\bm{x}_{0}$ and interate until $\left|\bm{x}_{i+1}-\bm{x}_{i}\right|/\left|\bm{x}_{i}\right|<\epsilon$, resulting in an $\bm{x}_{\text{\tiny{{WF}}}}$ that corresponds to the Wiener filter solution $\bm{s}_{\text{\tiny{{WF}}}}$. If convergence is not yet achieved, we re-initialise all parameters every 200 iterations and resume iterations, in order to facilitate convergence. An in depth review of the PCG algorithm is provided in Shewchuk (1994). Again, we stress the fact that finding a suitable preconditioner is the key factor when making use of the PCG method (e.g., Oh et al., 1999) and this consequently is the major stumbling block for state-of-the-art CMB data analysis.