# Wiener filter reloaded: fast signal reconstruction without preconditioning

Doogesh Kodi Ramanah, Guilhem Lavaux, Benjamin D. Wandelt

► To cite this version:

Doogesh Kodi Ramanah, Guilhem Lavaux, Benjamin D. Wandelt. Wiener filter reloaded: fast signal reconstruction without preconditioning. Monthly Notices of the Royal Astronomical Society, 2017, 468 (2), pp.1782-1793. 10.1093/mnras/stx527. hal-01554850

HAL Id: hal-01554850

https://hal.science/hal-01554850v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Wiener filter reloaded: fast signal reconstruction without preconditioning

Doogesh Kodi Ramanah,<fnref n="1" />,<fnref n="2" />* Guilhem Lavaux<fnref n="1" />,<fnref n="2" /> and Benjamin D. Wandelt<fnref n="1" />,<fnref n="2" />

Accepted 2017 February 27. Received 2017 February 23; in original form 2017 January 26

## 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 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 two and three to four times faster than the respective methods, for the specific problem considered.

**Key words:** methods: data analysis – methods: statistical – cosmic background radiation – cosmology: observations.

## 1 INTRODUCTION

In an era of precision cosmology, the data analysis of state-of-the-art galaxy redshift surveys and cosmic microwave background (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, Hoffman & Dekel 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. Eriksen et al. 2004; Jewell, Levin & Anderson 2004; O’Dwyer et al. 2004; Wandelt, Larson & Lakshminarayanan 2004; Larson et al. 2007; Smith, Zahn & Doré 2007; Elsner & Wandelt 2013; Bunn & Wandelt 2016).

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

$$d = s + n, \quad (1)$$

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

$$(S^{-1} + N^{-1})s_{\text{WF}} = N^{-1}d, \quad (2)$$

where $S$ and $N$ are the signal and noise covariance matrices, respectively, and $s_{\text{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 $S$ and $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

<footnote n="1">¹Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="2">²Sorbonne Universités, Institut Lagrange de Paris (ILP), 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: ramanah@iap.fr</footnote>
to the Wiener filter equation (e.g. Hirata et al. 2004; Komatsu, Spergel & Wandelt 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 (e.g. Eriksen et al. 2004; Wandelt 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, Spergel & Hinshaw 1999) since matrices are often extremely ill-conditioned in typical CMB problems (Eriksen et al. 2004). The inclusion of polarization 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 multilevel solver for Gaussian-constrained realizations 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.

This paper is organized 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 summarize 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 $t$, the so-called messenger field, with covariance $T$, where $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, that is, are sparse. As a result, although directly applying combinations like $(S + N)^{-1}$ may not be possible, we can always apply expressions like $(S + T)^{-1}$ and $(N + T)^{-1}$, irrespective of the basis chosen to render $S$ and $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 $t$, the modified $\chi^2$, where the posterior probability distribution of $s$ is proportional to $\exp(-\chi^2/2)$ is as follows:

$$ \chi_T^2 = (d - t)^{\dagger} \bar{N}^{-1} (d - t) + (t - s)^{\dagger} \mathbf{T}^{-1} (t - s) + s^{\dagger} \mathbf{S}^{-1} s, \quad (3) $$

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

$$ [\bar{N}^{-1} + (\lambda \mathbf{T})^{-1}] t = \bar{N}^{-1} d + (\lambda \mathbf{T})^{-1} s \quad (4) $$

$$ [\mathbf{S}^{-1} + (\lambda \mathbf{T})^{-1}] s = (\lambda \mathbf{T})^{-1} t, \quad (5) $$

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 initialize the vectors $s$ and $t$ with zeros, and choose an initial high value of $\lambda$. We first solve equation (4) for the messenger field $t$ in the basis defined by $N$, the noise covariance matrix. Then, we change to a basis where the signal covariance matrix $S$ has a sparse representation, such as Fourier space. Next, we solve for $s$ using equation (5) and the resulting $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. $s \to s_{WF}$, as $\lambda \to 1$.

If $\epsilon_i$ is the residual at the $i$-th step, then the following ($i+1$)-th step, the corresponding residual is

$$ \epsilon_{i+1} = [\mathbf{S}(\mathbf{S} + \lambda \mathbf{T})^{-1}] [\bar{N}(\bar{N} + \lambda \mathbf{T})^{-1}] \epsilon_i \quad (6) $$

$$ = [\mathbf{S}(\mathbf{S} + \lambda \mathbf{T})^{-1}] \{(\mathbf{N} - \mathbf{T})[\mathbf{N} + (\lambda - 1)\mathbf{T}]^{-1}\} \epsilon_i. \quad (7) $$

For the case of homogeneous noise, i.e. $N \propto \mathbf{1}$, the system is solved exactly in a single step. A key observation is that $|\epsilon_{i+1}| < |\epsilon_i|$ for all $i$ since the terms in brackets are less than unity, resulting in unconditional convergence of the signal reconstruction $s$ to the Wiener filter solution $s_{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 \le \eta \le 1$. The rationale behind the cooling scheme adopted in this work is laid out in Appendix A.
**Algorithm 1 Messenger algorithm**

1: **procedure** MESSENGER(d, N, S, N, L)

2: $s_0$ = zeros(N, N) ▷ Initialize s with zeros

3: $t_0 = d$ ▷ Initialize t via an initial guess

4: ▷ Compute the covariance of messenger field t

5: $\alpha$ = min(diag(N)) ▷ such that $T = \alpha \perp$

6: $\tilde{N} = N - T$ ▷ Compute the covariance $\tilde{N}$

7: **while** $\lambda = 10^4 \to 1$ **do**

8: **repeat**

9: ▷ Transform to Fourier space, $\mathcal{F}$

10: ▷ $\hat{s}_\ell = \mathcal{F}s_x = \left(\frac{L}{N}\right)^2 \sum_x \omega^{-x-k} s_x$

11: ▷ where $\omega = \exp(\frac{i2\pi\ell}{N})$

12: $\hat{s}_{i+1, \ell} = [\mathbf{S}^{-1} + \sigma(\lambda\mathbf{T})^{-1}]_{\ell}^{-1}$

13: $\mathcal{F}[(\lambda\mathbf{T})_x^{-1}t_{i,x}]_\ell$

14: ▷ $\sigma = \frac{N^2}{L^4}$ is a numerical factor due to $\mathcal{F}$

15: ▷ $N$ = no. of pixels, $L$ = angular extent

16: ▷ Transform to pixel space

17: $s_{i+1, x} = \mathcal{F}^{-1}(\hat{s}_{i+1, \ell})$

18: $t_{i+1, x} = [\tilde{\mathbf{N}}^{-1} + (\lambda \mathbf{T})^{-1}]_x^{-1}$

19: $\cdot [\tilde{\mathbf{N}}^{-1} \mathbf{d} + (\lambda \mathbf{T})^{-1} s_{i+1}]_x$

20: $i \leftarrow i + 1$

21: **until** $\|s_i - s_{i-1}\| / \|s_i\| < \epsilon$

22: $\lambda \leftarrow \lambda \times \eta$ ▷ Cooling scheme for $\lambda$

23: **end while**

24: $s \to s_{WF}$ ▷ as $\lambda = 1$

25: **return** $s_{WF}$

26: **end procedure**

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 are masked. Elsner & Wandelt (2012, 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 polarization 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 XXIII 2014; Planck Collaboration XXIV 2014). 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, Wandelt & Lavaux (2015) adopted a similar approach in their Bayesian hierarchical modelling of the CMB gravitational lensing. Alsing et al. (2016) 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, Heavens & Jaffe 2017). The works of Jasche & Lavaux (2015), Anderes et al. (2015) and Alsing et al. (2016) 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 optimized then becomes

$$ \chi_U^2 = (d-s)^\dagger \mathbf{N}^{-1} (d-s) + (s-u)^\dagger \mathbf{U}^{-1} (s-u) + u^\dagger \tilde{\mathbf{S}}^{-1} u, \quad (8) $$

where, analogous to the standard approach, $\mathbf{U} = v \perp$ with $v \equiv \min(\text{diag}(\mathbf{S}))$, and the covariance of the auxiliary field, $\tilde{\mathbf{S}} = \mathbf{S} - \mathbf{U}$. We derive the corresponding equations that must be satisfied by $s$ and $u$ at the minimum of $\chi_U^2$:

$$ (\mathbf{N}^{-1} + \mathbf{U}^{-1}) s = \mathbf{N}^{-1} d + \mathbf{U}^{-1} u \quad (9) $$

$$ (\mathbf{U}^{-1} + \tilde{\mathbf{S}}^{-1}) u = \mathbf{U}^{-1} s. \quad (10) $$

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

$$ \epsilon_{i+1} = [\mathbf{N}(\mathbf{N} + \mathbf{U})^{-1}][\tilde{\mathbf{S}}(\tilde{\mathbf{S}} + \mathbf{U})^{-1}]\epsilon_i \quad (11) $$

$$ = [\mathbf{N}(\mathbf{N} + \mathbf{U})^{-1}][\mathbf{S} - \mathbf{U}]\mathbf{S}^{-1}\epsilon_i. \quad (12) $$

This provides the basis for the following mechanism: we artificially truncate the spectrum $\mathbf{S}$ to some lower initial value of $\ell_{\text{iter}}$ that corresponds to a covariance $\mu$ and bring $\ell_{\text{iter}}$ slowly to $\ell_{\text{max}}$ corresponding to our final covariance $\nu$. So, essentially, we vary the covariance $\mathbf{U}$ via a cooling scheme to bring $\mu \to \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 $\tilde{\mathbf{S}}$ using the Heaviside function as $\tilde{\mathbf{S}} = \Theta(\mathbf{S} - \mathbf{U})$, as described quantitatively in Appendix B [cf. equations (B9) and (B10)]. For $\ell \lesssim \ell_{\text{iter}}$, the ratio given by equation (12) is always convergent. The algorithm consists of the following steps. We initialize the vectors $s$ and $u$ with zeros. We begin iterations with an initial value of $\ell_{\text{iter}}$ and correspondingly $\mu$, and iterate until $\mu \to \nu$ at $\ell_{\text{iter}} = \ell_{\text{max}}$, i.e. $s \to s_{\text{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 \to \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 \min(\text{diag}(\mathbf{N}))$, thereby
incorporating aspects of the standard messenger method into the dual messenger scheme. This leads to the following $\chi^2$:

$$ \chi_{\xi}^2 = (d - t)^{\dagger} \bar{\mathbf{N}}^{-1} (d - t) + (t - s)^{\dagger} \xi^{-1} (t - s) + s^{\dagger} \bar{\mathbf{S}}^{-1} s, \quad (13) $$

where $\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

$$ (\bar{\mathbf{N}}^{-1} + \xi^{-1}) t = \bar{\mathbf{N}}^{-1} d + \xi^{-1} s \quad (14) $$

$$ (\xi^{-1} + \bar{\mathbf{S}}^{-1}) s = \xi^{-1} t. \quad (15) $$

If $\alpha = 0$, we recover the usual dual messenger scheme, while setting $\mu = 0$ yields the standard messenger algorithm. The definition of $\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 $\xi$ by the factor $(1/\beta)$ and iterate until $\xi \to \alpha$, at which point $\mu = 0$ and we obtain the proper solution as desired.

### Algorithm 2 Dual messenger algorithm

1: **procedure** DUAL MESSENGER(d, N, S, N, L)

2: $s_0 = \text{zeros}(N, N)$ ▷ Initialize $s$ with zeros

3: $t_0 = d$ ▷ Initialize $t$ via an initial guess

4: ▷ Compute the covariance of auxiliary field $t$

5: $\alpha = \min(\text{diag}(\mathbf{N}))$ ▷ such that $\mathbf{T} = \alpha\mathbb{1}$

6: $\bar{\mathbf{N}} = \mathbf{N} - \mathbf{T}$ ▷ Compute the covariance $\bar{\mathbf{N}}$

7: **while** $\xi = (\alpha + \mu) \to \alpha$ **do**

8: $\mathbf{U} = (\sigma \mu)\mathbb{1}$ ▷ Compute covariance $\mathbf{U}$

9: ▷ As in Algorithm 1, factor of $\sigma$ due to $\mathcal{F}$

10: $\bar{\mathbf{S}} = \Theta(\mathbf{S} - \mathbf{U})$ ▷ Compute covariance $\bar{\mathbf{S}}$

11: **repeat**

12: ▷ Moving to Fourier space, $\mathcal{F}$

13: $\hat{t}_{i+1,\ell} = [(\bar{\mathbf{S}}^{-1} + \sigma\xi^{-1})^{-1}\sigma\xi^{-1}]_{\ell} \mathcal{F}(t_{i,x})$

14: ▷ Transform to pixel space

15: $s_{i+1,x} = \mathcal{F}^{-1}(\hat{t}_{i+1,\ell})$

16: $t_{i+1,x} = (\mathbf{N}^{-1} + \xi^{-1})_x^{-1}$

17: $\cdot (\mathbf{N}^{-1} d + \xi^{-1} s_{i+1})_x$

18: $i \leftarrow i + 1$

19: **until** $\|s_i - s_{i-1}\| / \|s_i\| < \epsilon$

20: $\xi \leftarrow \xi \times \beta$ ▷ Cooling scheme for $\xi$

21: $\mu \leftarrow (\xi - \alpha)/\sigma$ ▷ Compute resulting $\mu$

22: **end while**

23: $s \to s_{\text{WF}}$ ▷ as $\xi = \alpha, \mu = 0$

24: **return** $s_{\text{WF}}$

25: **end procedure**

## 3 APPLICATION TO CMB

We generated artificial CMB data, in a realistic scenario, by drawing Gaussian random fields on a 2D flat sky with grid resolution, $N_{\text{pixels}} = 512^2$, and angular extent, $L = 10.0$ deg. This true simulated map is contaminated by Gaussian white noise with covariance matrix, $\mathbf{N} = (64.0 \mu K^2)\mathbb{1}$, via a linear model as given in equation (1). We made use of CAMB<fnref n="1" /> (Lewis, Challinor & Lasenby 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 XIII 2016). A portion of the flat sky is masked by imposing a central square patch of angular extent 5.0 deg with a noise covariance of $10^6$ 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, $\|s_{i+1} - s_i\| / \|s_i\| < \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 Fig. 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 Fig. 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 Fig. 2, 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 per cent residuals. The dual messenger and PCG reconstructions yield around 4 and 9 per cent 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 Fig. 3, 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 Fig. 3. Standard Wiener filtering algorithms usually encounter difficulties in dealing with masked regions having infinite noise. By masking 25 per cent of the simulated map with noise of the order of $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 polarization 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 polarization components of the signal covariance. The Wiener filtering of polarized CMB data with more complex noise models

<footnote n="1">¹ http://camb.info</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The simulated CMB map and the reconstruction via the dual messenger algorithm. The left-hand panel shows the simulated map, which is subsequently contaminated by white noise, with the central square patch masked by extremely high-noise covariance. The right-hand panel illustrates the corresponding reconstructed map obtained via the dual messenger technique. The messenger and PCG reconstructions are not shown, as they look similar to the dual messenger reconstructed map. The residual maps displayed in Fig. 2, however, help to discern the differences between the various reconstructions.</figcaption></figure>

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.

Fig. 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 \to 1$ in a three-step procedure, as depicted in Fig. 4 by 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_{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 Fig. 4 by dashed blue lines.

Fig. 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 s, whereas the messenger and PCG methods have corresponding wall-clock times of roughly 435 and 663 s, 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, $O(N_{pixels} \log N_{pixels})$, and two scalar multiplications corresponding to algebraic operations of $O(N_{pixels})$, per iteration. In terms of memory requirements, two vectors of size $N_{pixels}$ must be temporarily stored in memory. In comparison, the PCG method requires three Fourier transforms and 10 scalar multiplications per iteration, and temporary storage of eight vectors of dimension $N_{pixels}$ in memory.

We also performed a consistency check by verifying the variation of the residual error given by $\|Ax - y\| / \|y\|$. This is usually adopted as a convergence criterion for PCG computations. The variation of this residual error, analogous to Fig. 4, is illustrated in Fig. 5. The corresponding residual errors obtained via the different algorithms drop below their respective Cauchy convergence thresholds implemented. The oscillatory behaviour of the PCG solution, displayed in Figs 4 and 5, is mainly due to the reinitialization step after every 200 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 Fig. 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 (A5) and (B14)], 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.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The reference map and residual maps from the three algorithms. The top left-hand panel depicts the reference map computed using the PCG scheme with $\epsilon = 10^{-9}$, while the other panels illustrate the residual maps yielded by the three methods, generated by computing the difference between the reference map and the corresponding reconstructed maps over the full sky. The messenger approach produces the least amount of residuals, around 3 per cent, while the dual messenger and PCG schemes result in approximately 4 and 9 per cent residuals, respectively. For the messenger reconstructions, the residuals lie mostly on the edges of the mask.</figcaption></figure>

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 Fig. 6, with the left-hand and right-hand 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_M^2$ 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 fig. 2 of Elsner & Wandelt (2012). The $\chi_{DM}^2$ of the dual messenger algorithm has a similar behaviour. The $\chi^2$ in both cases finally attains the $\chi_{ref}^2$ of the reference PCG method with $\epsilon = 10^{-9}$, with $\Delta\chi_M^2/\chi_{ref}^2 = 9.2 \times 10^{-6}$ and $\Delta\chi_{DM}^2/\chi_{ref}^2 = 9.0 \times 10^{-6}$, where $\chi_{ref}^2 = 4.2 \times 10^4$.

Another important convergence diagnostic is the variation of the relative error, $C_\ell(s_{WF} - s_{WF}^{\text{ref}})/C_\ell(s_{WF}^{\text{ref}})$, computed over the full sky, as a function of scale, $\ell$, depicted in Figs 7, 8 and 9 for the messenger, dual messenger and PCG algorithms, respectively. The relative errors on the small scales are of the order of $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. Fig. 2). 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 Fig. 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,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Reconstructed power spectra computed using the different algorithms. The dashed line indicates the input angular power spectrum, $C_{\ell}^{\text{input}}$, from which the CMB signals are drawn. The power spectra recovered by the three algorithms are all in good agreement with the reference power spectrum computed using PCG method with $\epsilon = 10^{-9}$ on all scales. The deviations on the small scales from the input power spectrum are due to the characteristic feature of a Wiener-filtered signal, where the power on small scales, in the low signal-to-noise regime, is suppressed.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variation of the residual error, given by the Cauchy criterion, with number of iterations for the different methods. The Cauchy convergence criteria $\epsilon_i$ imposed for the different regimes of the cooling schemes for $\lambda$ and $\mu$ are also displayed, along with the corresponding thresholds $\epsilon$ used for the two PCG methods, in dashed lines. The vertical green line denotes the convergence point of the PCG scheme, with its convergence behaviour already represented by the reference PCG method. The dual messenger algorithm requires the smallest number of iterations to converge to the final solution than the other two methods.</figcaption></figure>

and therefore computation time, by around 25 per cent 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

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Variation of the residual error, given by $\|Ax - y\| / \|y\|$, with number of iterations for the different algorithms. The corresponding residual errors for all methods drop below the convergence thresholds adopted, demonstrating the consistency of our computations dictated by the Cauchy criteria (cf. Fig. 4).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Variation of $\chi^2$ with number of iterations for the two messenger methods. In the left-hand panel, we find that the $\chi_M^2$ of the messenger algorithm drops rapidly with the cooling scheme for $\lambda$ adopted, as expected from the earlier discussion. In the right-hand panel, the $\chi_{DM}^2$ of its dual counterpart displays a similar behaviour. In both cases, the final solution matches the $\chi_{ref}^2$ of the PCG method with $\epsilon = 10^{-9}$.</figcaption></figure>

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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Convergence rate by frequency bin for the messenger method. This illustrates the relative error as a function of scale, $\ell$. Each line in the figure corresponds to the relative error for a specific value of the scalar parameter $\lambda$. The messenger algorithm converges smoothly and in nearly uniform fashion on small and intermediate scales, with the relative error dropping till below $10^{-6}$. However, for larger scales, the relative error is reduced by lower extent, but stays below $10^{-3}$ for the largest scales. For $\lambda = 1$, the behaviour is similar to that displayed by PCG (cf. Fig. 9).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Convergence rate by frequency bin for the dual messenger method. Same as Fig. 7, except that here each line in the figure corresponds to a given value of $\xi$, displaying the hierarchical nature of this scheme. The final relative error is sufficiently low on all scales. A quantitative explanation of the convergence behaviour is given in Appendix B.</figcaption></figure>

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 optimized 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 threshold-

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Convergence rate by frequency bin for the PCG method. Same as Fig. 7. The relative error as a function of scale for every 100th iteration is illustrated. The PCG method has a lower relative error than the dual messenger algorithm on the small scales, but the behaviour on the largest scales remains inferior to that of both messenger methods. The relative spacing of the lines indicates that convergence slows down as the PCG method progresses in iterations, and this is in agreement with its variation of the residual error with number of iterations, as displayed by the green line in Fig. 4.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Execution times required for convergence as a function of convergence criterion $\epsilon$ for the different methods. This demonstrates the efficient performance of the two messenger algorithms with respect to the PCG method.</figcaption></figure>

ing 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_{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 generalization of both the messenger and dual messenger algorithms in a combined approach that may
further optimize the performance and efficiency compared to traditional techniques of solving the Wiener filter problem.

## ACKNOWLEDGEMENTS

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

## REFERENCES

Alsing J., Heavens A., Jaffe A. H., Kiessling A., Wandelt B., Hoffmann T., 2016, MNRAS, 455, 4452

Alsing J., Heavens A. F., Jaffe A. H., 2017, MNRAS, 466, 3272

Anderes E., Wandelt B. D., Lavaux G., 2015, ApJ, 808, 152

Bunn E. F., Wandelt B., 2016, preprint (arXiv:1610.03345)

Elsner F., Wandelt B. D., 2012, preprint (arXiv:1211.0585)

Elsner F., Wandelt B. D., 2013, A&A, 549, A111

Erdoğdu P. et al., 2004, MNRAS, 352, 939

Erdoğdu P. et al., 2006, MNRAS, 373, 45

Eriksen H. K. et al., 2004, ApJS, 155, 227

Eriksen H. K., Dickinson C., Jewell J. B., Banday A. J., Górski K. M., Lawrence C. R., 2008, ApJ, 672, L87

Gauss C. F., 1809, Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium Auctore Carolo Friderico Gauss. Frid. Perthes et IH Besser, Hamburg

Golub G. H., Van Loan C. F., 1996, Matrix Computations. Johns Hopkins Univ. Press, Baltimore, MD, p. 374

Hirata C. M., Padmanabhan N., Seljak U., Schlegel D., Brinkmann J., 2004, Phys. Rev. D, 70, 103501

Jasche J., Lavaux G., 2015, MNRAS, 447, 1204

Jasche J., Kitaura F. S., Wandelt B. D., Enßlin T. A., 2010, MNRAS, 406, 60

Jaynes E. T., 1957, Phys. Rev., 106, 620

Jewell J., Levin S., Anderson C. H., 2004, ApJ, 609, 1

Kalman R. E., 1960, J. Basic Eng., 82, 35

Kitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497

Kitaura F. S., Jasche J., Li C., Enßlin T. A., Metcalf R. B., Wandelt B. D., Lemson G., White S. D. M., 2009, MNRAS, 400, 183

Komatsu E., Spergel D. N., Wandelt B. D., 2005, ApJ, 634, 14

Larson D. L., Eriksen H. K., Wandelt B. D., Górski K. M., Huey G., Jewell J. B., O'Dwyer I. J., 2007, ApJ, 656, 653

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Mangilli A., Verde L., 2009, Phys. Rev. D, 80, 123007

Mangilli A., Wandelt B., Elsner F., Liguori M., 2013, A&A, 555, A82

O'Dwyer I. J., Eriksen H. K., Wandelt B. D., Jewell J. B. et al., 2004, ApJ, 617, L99

Oh S. P., Spergel D. N., Hinshaw G., 1999, ApJ, 510, 551

Planck Collaboration XXIII, 2014, A&A, 571, A23

Planck Collaboration XXIV, 2014, A&A, 571, A24

Planck Collaboration XIII, 2016, A&A, 594, A13

Seljebotn D. S., Mardal K.-A., Jewell J. B., Eriksen H. K., Bull P., 2014, ApJS, 210, 24

Shewchuk J. R., 1994, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. Department of Computer Science, Carnegie-Mellon University, Pittsburgh, PA

Smith K. M., Zahn O., Doré O., 2007, Phys. Rev. D, 76, 043510

Wandelt B. D., Larson D. L., Lakshminarayanan A., 2004, Phys. Rev. D, 70, 083511

Wiener N., 1949, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, Vol. 2. MIT Press, Cambridge, MA

Zaroubi S., 2002, MNRAS, 331, 901

Zaroubi S., Hoffman Y., Dekel A., 1999, ApJ, 520, 413

## APPENDIX A: COOLING SCHEME FOR MESSENGER ALGORITHM

In this appendix, 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 \mathbf{1}$, where $\alpha = \min(\text{diag}(\mathbf{N}))$, the $\chi_T^2$ for the messenger scheme from equation (3) can be written as

$$ \chi_T^2 = (\mathbf{d} - \mathbf{s})^\dagger [\mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} (\mathbf{d} - \mathbf{s}) + s^\dagger \mathbf{S}^{-1} s. \quad (A1) $$

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

$$ \begin{aligned} s_{\text{WF}}(\lambda) &= \{[\mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} + \mathbf{S}^{-1}\}^{-1} [\mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} \mathbf{d} \\ &= \mathbf{S} [\mathbf{S} + \mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} \mathbf{d}, \end{aligned} \quad (A2) $$

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

$$ \begin{aligned} s_{\text{WF}}(\lambda) - s_{\text{WF}}(\lambda') &= \mathbf{S} \{[\mathbf{S} + \mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} - [\mathbf{S} + \mathbf{N} + (\lambda' - 1)\alpha \mathbf{1}]^{-1}\} \mathbf{d} \\ &= \mathbf{S}[\mathbf{S} + \mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} \{\mathbf{S} + \mathbf{N} + (\lambda' - 1)\alpha \mathbf{1} - [\mathbf{S} + \mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]\} [\mathbf{S} + \mathbf{N} + (\lambda' - 1)\alpha \mathbf{1}]^{-1} \mathbf{d} \\ &= \alpha(\lambda' - \lambda)\mathbf{S}[\mathbf{S} + \mathbf{N} + (\lambda - 1)\alpha \mathbf{1}]^{-1} [\mathbf{S} + \mathbf{N} + (\lambda' - 1)\alpha \mathbf{1}]^{-1} \mathbf{d}. \end{aligned} \quad (A3) $$

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

$$ \|s_{\text{WF}}(\lambda)\|_{\ell} = C_{\ell}(C_{\ell} + \lambda\alpha)^{-1} \|d\|_{\ell}. \quad (A4) $$

Using equations (A3) and (A4), the relative error can then be computed as follows:

$$ \frac{\|s_{\text{WF}}(\lambda) - s_{\text{WF}}(\lambda')\|_{\ell}}{\|s_{\text{WF}}(\lambda)\|_{\ell}} \le \left[ |\lambda' - \lambda|\alpha \left( \frac{C_{\ell}}{C_{\ell} + \lambda\alpha} \right) \left( \frac{1}{C_{\ell} + \lambda'\alpha} \right) \|d\|_{\ell} \right] \left[ \frac{C_{\ell}}{C_{\ell} + \lambda\alpha} \|d\|_{\ell} \right]^{-1} \le \frac{\alpha |\lambda' - \lambda|}{C_{\ell} + \lambda'\alpha}, \quad (A5) $$
where, for an arbitrary matrix **M**, ||**M**||_ℓ is the subset of ||**M**|| over subspace ℓ. From equation (6), we have the fractional error reduction, after one iteration, given by

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

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

$$\epsilon \le 1 - \frac{C_{\ell}}{C_{\ell} + \lambda\alpha} = \frac{\lambda\alpha}{C_{\ell} + \lambda\alpha}. \qquad (A7)$$

To favour convergence, from equations (A5) and (A7), we require the two conditions:

$$\frac{\alpha |\lambda' - \lambda|}{C_{\ell} + \lambda' \alpha} \lesssim \begin{cases} \frac{\lambda \alpha}{C_{\ell} + \lambda \alpha} \\ \frac{\lambda' \alpha}{C_{\ell} + \lambda' \alpha} \end{cases} \qquad (A8)$$

so that the error reduction in one iteration matches the change in solution arising from two consecutive values of λ. This results in the following two constraints: |λ' - λ| ≤ λ' and |λ' - λ| ≤ λ. The former is automatically satisfied, while the latter leads to the following bounds: λ'/2 ≤ λ ≤ λ'. Therefore, for the cooling scheme, we can write λ = ηλ', where 1/2 ≤ η ≤ 1, to improve convergence on all scales. This is the motivation behind our choice of η = 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 (A2) for this scheme. We begin by writing down equations (9) and (10), showing the dependence on the power spectrum truncation μ:

$$(\mathbf{N}^{-1} + \mathbf{U}_{\mu}^{-1}) s_{\text{WF}}(\mu) = \mathbf{N}^{-1} d + \mathbf{U}_{\mu}^{-1} \mathbf{u}_{\mu} \qquad (B1)$$

$$(\mathbf{U}_{\mu}^{-1} + \bar{\mathbf{S}}_{\mu}^{-1}) \mathbf{u}_{\mu} = \mathbf{U}_{\mu}^{-1} s_{\text{WF}}(\mu). \qquad (B2)$$

Solving for $s_{\text{WF}}(\mu)$ via the following steps:

$$\begin{aligned} \mathbf{N}^{-1} d &= [(\mathbf{N}^{-1} + \mathbf{U}_{\mu}^{-1}) - \mathbf{U}_{\mu}^{-1}(\mathbf{U}_{\mu}^{-1} + \bar{\mathbf{S}}_{\mu}^{-1})^{-1}\mathbf{U}_{\mu}^{-1}] s_{\text{WF}}(\mu) = [(\mathbf{N}^{-1} + \mathbf{U}_{\mu}^{-1})(\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu}) - \bar{\mathbf{S}}_{\mu}\mathbf{U}_{\mu}^{-1}] (\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu})^{-1} s_{\text{WF}}(\mu) \\ &= [\mathbf{N}^{-1}(\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu}) + 1] (\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu})^{-1} s_{\text{WF}}(\mu) = [\mathbf{N}^{-1} + (\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu})^{-1}] s_{\text{WF}}(\mu), \end{aligned} \qquad (B3)$$

so that finally, we have

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

In the limit $\mu \to \nu$, where $\nu \equiv \min(\text{diag}(\mathbf{S}))$, $\bar{\mathbf{S}} + \mathbf{U} \to \mathbf{S}$, reducing equation (B4) to the usual Wiener filter equation (2), implying consistency. When we change the truncation of the spectrum, $\mu \to \tilde{\mu}$, for $\tilde{\mu} < \mu$, we obtain a new solution $s_{\text{WF}}(\tilde{\mu})$. The counterpart of equation (A5) is then

$$\begin{aligned} s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu) &= (\bar{\mathbf{S}}_{\tilde{\mu}} + \mathbf{U}_{\tilde{\mu}}) [(\bar{\mathbf{S}}_{\tilde{\mu}} + \mathbf{U}_{\tilde{\mu}}) + \mathbf{N}]^{-1} d - (\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu}) [(\bar{\mathbf{S}}_{\mu} + \mathbf{U}_{\mu}) + \mathbf{N}]^{-1} d \\ &= [(\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}] d. \end{aligned} \qquad (B5)$$

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

$$\bar{\mathbf{S}}_{\tilde{\mu}} = \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu}, \mu} \mathbf{1}, \qquad (B6)$$

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:

$$\bar{\mathbf{S}}_{\mu} = \Theta(\mathbf{S} - \mu) \qquad (B7)$$

$$\bar{\mathbf{S}}_{\tilde{\mu}} = \Theta(\mathbf{S} - \tilde{\mu}), \qquad (B8)$$

where, for a matrix **M** = **P**Λ**P**⁻¹, after applying a basis transformation, with Λ being diagonal,

$$\Theta(\mathbf{M}) = \mathbf{P}\Theta(\Lambda)\mathbf{P}^{-1}, \qquad (B9)$$

and

$$\Theta(\Lambda)_{ii} = \begin{cases} 0, & \Lambda_{ii} \le 0 \\ \Lambda_{ii}, & \Lambda_{ii} > 0. \end{cases} \qquad (B10)$$
Using equation (B6) in equation (B5) results in

$$
\begin{align*}
s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu) &= [(\bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})^{-1} - (\bar{\mathbf{S}}_{\mu} - \mu)(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}] \mathbf{d} \\
&= [\bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu} - (\bar{\mathbf{S}}_{\mu} + \mu)(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})] (\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu + \Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)^{-1} \mathbf{d} \\
&= [\bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \bar{\mathbf{S}}_{\mu} - \mu - (\bar{\mathbf{S}}_{\mu} + \mu) + (\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)] (\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})^{-1} \mathbf{d} \\
&= [(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu) - \bar{\mathbf{S}}_{\mu} - \mu] (\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})^{-1} \mathbf{d} \\
&= \mathbf{N}(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \Delta_{\tilde{\mu},\mu} + \tilde{\mu})^{-1} \mathbf{d}. \tag{B11}
\end{align*}
$$

Using equation (B6) and substituting the following form of equation (B4),

$$
s_{\text{WF}}(\tilde{\mu}) = \bar{S}_{\tilde{\mu}} [\mathbf{N} + (\bar{S}_{\tilde{\mu}} + \tilde{\mu})^{-1}]^{-1} \mathbf{d}, \quad (\text{B12})
$$

in equation (B11) leads to

$$
s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu) = \mathbf{N}(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)(\bar{\mathbf{S}}_{\tilde{\mu}} + \tilde{\mu})^{-1}s_{\text{WF}}(\tilde{\mu}), \quad \text{(B13)}
$$

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

$$
\frac{\|s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu)\|}{\|s_{\text{WF}}(\tilde{\mu})\|} \le \frac{\|\mathbf{N}(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}(\Delta_{\tilde{\mu},\mu} + \tilde{\mu} - \mu)(\bar{\mathbf{S}}_{\tilde{\mu}} + \tilde{\mu})^{-1}\|}{\|\bar{\mathbf{S}}_{\tilde{\mu}} + \tilde{\mu}\|} \quad \text{(B14)}
$$

Now, for $\mu \to \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 (B7) and (B8) in equation (B6),

$$
\Theta(\mathbf{S} - \tilde{\mu}) = \Theta(\mathbf{S} - \mu) + \Delta_{\tilde{\mu},\mu}. \quad (B15)
$$

The three regimes are:

(i) $\mathbf{S}_{\ell} > \mu \mathbf{1}$,

where $\Theta(\mathbf{S} - \tilde{\mu}) = \mathbf{S}_{\ell} - \tilde{\mu} \mathbf{1}$, and $\Theta(\mathbf{S} - \mu) = \mathbf{S}_{\ell} - \mu \mathbf{1}$, so that $\Delta_{\tilde{\mu},\mu} = \mu - \tilde{\mu}$, and this causes the second term in equation (B14) 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.

(ii) $\tilde{\mu} \mathbf{1} < \mathbf{S}_{\ell} < \mu \mathbf{1}$,

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

$$
\frac{\|s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu)\|}{\|s_{\text{WF}}(\tilde{\mu})\|} \le \frac{\|\mathbf{N}(\mathbf{N} + \bar{\mathbf{S}}_{\mu} + \mu)^{-1}\|}{\|\mathbf{S}_{\ell} - \mu \mathbf{1}\|} \|(\bar{\mathbf{S}}_{\tilde{\mu}} + \tilde{\mu})^{-1}\|_{\tilde{\mu},\mu}. \quad (\text{B16})
$$

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

$$
\frac{\|s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu)\|}{\|s_{\text{WF}}(\tilde{\mu})\|} \le \alpha' \left| \frac{C_{\ell} - \mu}{C_{\ell}} \right|_{\tilde{\mu},\mu} = \alpha' \left| 1 - \frac{\mu}{C_{\ell}} \right|_{\tilde{\mu},\mu}, \quad (B17)
$$

since $\bar{\mathbf{S}}_{\tilde{\mu}} + \tilde{\mu} \sim C_{\ell}$. Also, $\alpha' \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_{\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.

(iii) $\mathbf{S}_{\ell} < \tilde{\mu} \mathbf{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_{\ell}$ yields

$$
\frac{\|s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu)\|}{\|s_{\text{WF}}(\tilde{\mu})\|} \le \alpha' \left| \frac{\tilde{\mu} - \mu}{C_{\ell}} \right|_{\tilde{\mu},\mu}. \quad (B18)
$$

So, the overall convergence behaviour can be quantitatively described by

$$
\frac{\|s_{\text{WF}}(\tilde{\mu}) - s_{\text{WF}}(\mu)\|}{\|s_{\text{WF}}(\tilde{\mu})\|} \le 0 + \alpha' \left| 1 - \frac{\mu}{C_{\ell}} \right|_{\tilde{\mu},\mu} + \alpha' \left| \frac{\tilde{\mu} - \mu}{C_{\ell}} \right|_{\tilde{\mu},\mu}, \quad (B19)
$$

and this leads to the following interpretation: for higher values of $\mu$ on large scales, the third term in equation (B19) dominates since $|\tilde{\mu} - \mu|$ 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_\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 PCG METHOD

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

$$A\mathbf{x} = \mathbf{y}, \tag{C1}$$

where $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 $M$, which is the approximate inverse of $A$, i.e. $M \approx A^{-1}$, as follows:

$$M A \mathbf{x} = M \mathbf{y}. \tag{C2}$$

The Wiener filter equation (2) can be rewritten as

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

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

$$M_{\ell\ell}^{-1} = A_{\ell\ell} = \left(\sum_k n_k\right) C_{\ell} \frac{1}{L^4}, \tag{C4}$$

where $n_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 $\mathbf{x}_0$ and iterate until $|\mathbf{x}_{i+1} - \mathbf{x}_i| / |\mathbf{x}_i| < \epsilon$, resulting in an $\mathbf{x}_{\text{WF}}$ that corresponds to the Wiener filter solution $\mathbf{s}_{\text{WF}}$. If convergence is not yet achieved, we reinitialize 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 is consequently the major stumbling block for state-of-the-art CMB data analysis.

This paper has been typeset from a TeX/LATEX file prepared by the author.