Preconditioner-free Wiener filtering with a dense noise matrix

TL;DR

This paper introduces a new iterative method for Wiener filtering that efficiently handles dense noise covariance matrices by decomposing them into sparse components, enabling applications to complex Cosmic Microwave Background data.

Contribution

It extends the Elsner & Wandelt (2013) method to dense noise matrices using multiple messenger fields, allowing for efficient Wiener filtering in more realistic scenarios.

Findings

Successfully reproduces Wiener filter solutions for test problems.

Applies the method to a simulated CMB map with complex noise properties.

Addresses challenges in filtering data from advanced ground-based CMB experiments.

Abstract

This work extends the Elsner & Wandelt (2013) iterative method for efficient, preconditioner-free Wiener filtering to cases in which the noise covariance matrix is dense, but can be decomposed into a sum whose parts are sparse in convenient bases. The new method, which uses multiple messenger fields, reproduces Wiener filter solutions for test problems, and we apply it to a case beyond the reach of the Elsner & Wandelt (2013) method. We compute the Wiener filter solution for a simulated Cosmic Microwave Background map that contains spatially-varying, uncorrelated noise, isotropic $1/f$ noise, and large-scale horizontal stripes (like those caused by the atmospheric noise). We discuss simple extensions that can filter contaminated modes or inverse-noise filter the data. These techniques help to address complications in the noise properties of maps from current and future generations of ground-based Microwave Background experiments, like Advanced ACTPol, Simons Observatory, and CMB-S4.

Equations

$$d=s_{\mathrm{true}}+n.$$

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

$$s=({\mathbf{S}}^{-1}+{\mathbf{N}}^{-1}{)}^{-1}{\mathbf{N}}^{-1}d.$$

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

$$t$$

$$s$$

$$t_1$$

$$d_0$$

$$t_0$$

$$s$$

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

$$d_{\mathrm{rec}}(s)=(\mathbf{S}+\mathbf{N}){\mathbf{S}}^{-1}s=(\mathbf{I}+\mathbf{N}{\mathbf{S}}^{-1})s.$$

$$\Delta d=d_{\mathrm{rec}}-d,$$

$$d=s_{\mathrm{CMB}}+n_0+n_1.$$

$$⟨n_{0,lm}{n}_{0,l^{\prime }{m}^{\prime }}^{\ast }⟩=N_{0,lm}{\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$${\mathbf{N}}_i=\mathbf{F\Lambda }{\mathbf{F}}^{†},$$

$${\bar{\mathbf{N}}}_i^{-1}$$

$${\mathbf{N}}^{-1}d\approx {\alpha }^{-1}\left[\alpha \mathbf{I}(\alpha \mathbf{I}+\mathbf{N}{)}^{-1}d\right],$$

$$-2\log P(s,d)=s^{†}{\mathbf{S}}^{-1}s+(d-s{)}^{†}{\mathbf{N}}^{-1}(d-s)+\text{mbox}const.$$

$$-2\log P(s|d)=(s-{\mu }_s{)}^{†}({\mathbf{S}}^{-1}+{\mathbf{N}}^{-1})(s-{\mu }_s)+\text{mbox}const.$$

$$-2\log P_1(s,d,t)$$

$$\int dt{P}_1(s,d,t)=P(s,d).$$

$$t$$

$$s$$

$$\mathbf{S}$$

$$-2\log P_2(s,t_0,d_0,t_1,d)$$

$${\mathbf{T}}_0$$

$${\mathbf{T}}_1$$

$$t_0=[{\bar{\mathbf{N}}}_0^{-1}+{\mathbf{T}}_0^{-1}-{\mathbf{N}}_0^{-1}({\bar{\mathbf{N}}}_0^{-1}+{\mathbf{T}}_1^{-1}){\mathbf{N}}_0^{-1}{]}^{-1}$$

$$\left({\bar{\mathbf{N}}}_0^{-1}({\bar{\mathbf{N}}}_0^{-1}+{\mathbf{T}}_0^{-1}{)}^{-1}{\mathbf{T}}_1^{-1}{t}_1+{\mathbf{T}}_0^{-1}s\right),$$

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

Full text

1 Introduction

The Wiener filter (Wiener, 1949) is a general method to estimate a stochastic field based on noisy or incomplete data. It has long been used in astrophysics and cosmology to reconstruct sparse data (e.g. Rybicki & Press, 1992), denoise maps of the Cosmic Microwave Background (CMB, e.g. Bunn et al., 1994) and large scale structure (e.g. Zaroubi et al., 1995), estimate power spectra (e.g. Oh et al., 1999), and combine lensing and direct probes of the large scale structure (e.g. Simon, 2013; Alsing et al., 2016), among other uses. It closely relates to Bayesian methods for map and power spectrum estimation (e.g. Gibbs sampling, Wandelt et al., 2004; Eriksen et al., 2004).

The Wiener filter examines a data vector ($d$) that is a linear combination of signal ($s_{\rm true}$) and stochastic noise ($n$):

[TABLE]

If the signal covariance is $\mathbf{S}$ and the noise covariance is $\mathbf{N}$, the Wiener filter estimate for the signal is

[TABLE]

or equivalently

[TABLE]

For a particular realization of the data, the Wiener filter solution minimizes:

[TABLE]

maximizing the multivariate Gaussian probability $\propto\exp(-\chi^{2}/2)$. No other linear solution has reconstruction errors $\epsilon=s-s_{\rm true}$ with smaller mean-square deviations: the Wiener filter minimizes $\langle\epsilon^{\dagger}\epsilon\rangle$, averaged over all signal and noise realizations.

For large data sets, the solution is difficult unless both $\mathbf{S}$ and $\mathbf{N}$ matrices are sparse (making them easier to invert) and sparse in the same basis (making the sums in equations 2 and 3 easier to tabulate and invert). In cosmological contexts, the signal covariance is often isotropic, and so is diagonal in harmonic space. The noise covariance, by contrast, is often linked to the survey strategy, and so may be dominantly diagonal in real space. This mismatch has meant that Wiener filter solutions have traditionally required sophisticated linear system solvers that avoid the inversions of dense matrices, using for example iterative conjugate gradient methods (Hirata et al., 2004; Smith et al., 2007), and these require careful attention to their preconditioners.

Elsner & Wandelt (2013) circumvented this difficulty by devising an iterative method to solve the Wiener filter equation without a preconditioner in cases where the signal and noise are sparse in different bases. To do so they introduce an auxiliary messenger field, which splits off a homogeneous component of the noise covariance. Such a component is diagonal in every orthogonal basis. This general formalism provides new tools to tackle problems with isotropic signal and spatially varying but uncorrelated noise, and also closely relates to a method for Gibbs sampling. Already various authors have employed these methods in the context of CMB polarization (Elsner & Wandelt, 2013), CMB gravitational lensing (Anderes et al., 2015), large-scale structure (Jasche & Lavaux, 2015; Lavaux & Jasche, 2016), and cosmic shear (Alsing et al., 2016). In all of these cases the signal covariance is (block-) diagonal in harmonic space and the noise covariance is (block-) diagonal in pixel space. These methods are under rapid development, and recently Kodi Ramanah et al. (2017) showed that additionally splitting off a homogeneous portion of the signal covariance can help to speed the convergence.

In this work we extend this type of solution to a new class of problems. We consider cases in which the noise covariance is dense, but is a sum of pieces that can each be represented sparsely in convenient bases. Our approach achieves this by adding messenger fields for each additional term required for the covariance.

This paper is organized so that in section 2 we review the Elsner & Wandelt (2013) method and introduce our extension. In section 3 we apply our method to simulations with progressively more complex noise. In section 4 we discuss the implications. An appendix provides more details of our approach.

2 Methods

Elsner & Wandelt (2013) (hereafter EW) showed that the Wiener filter equation could be recast as a system of equations with the same solution. For a matrix $\mathbf{T}$ (which is arbitrary except that $\mathbf{\bar{N}=N-T}$ should be positive definite), the solution for the Wiener filter also solves

[TABLE]

for a unique $t$. (The value of $t$ depends on our choice for $\mathbf{T}$.) Furthermore, iterating equations (5)–(6) will exponentially converge to the Wiener solution. Here we solve for two vectors ($s,t$) instead of one, which makes the problem look harder until we notice that we can make specific, convenient choices for $\mathbf{T}$. If $\mathbf{T}$ is proportional to the identity matrix ($\mathbf{T}=\tau_{\rm EW}\mathbf{I}$, where $\tau_{\rm EW}$ is a scalar constant) then it is sparse in every orthonormal basis. Thus if $\mathbf{S}$ and $\mathbf{N}$ can separately be written in distinct sparse bases, a combination like $(\mathbf{S}^{-1}+\mathbf{T}^{-1})$ can be written in a sparse way and easily inverted. That allows the equations to be evaluated directly without recourse to linear system solvers. This is especially quick if the signal and data vectors can be transformed to the convenient bases via a fast transform (e.g. between real and harmonic space). A “cooling schedule” that artificially scales up the $\mathbf{T}$ covariance at the start, then gradually returns it to the proper value, can help speed the convergence.

Although the EW equations (5)–(6) are always valid, it only makes sense to apply them in cases when both the signal and noise have convenient sparse representations that allow the above shortcut. However, for many problems of interest, the noise covariance is dense for all convenient bases. One example problem is the optimal filtering of maps from ground-based CMB experiments, where a significant amount of the noise comes from the atmosphere. That portion of the noise covariance is strongly correlated in real space between different parts of the map, while the portion of the noise covariance due to detector noise and foreground masking are sparse only in real space. Thus the noise covariance is dense both in real and harmonic space and is by itself difficult to invert.

Here we show that we can make progress in this case if we can decompose the noise covariance into a sum of simpler terms. If each term in the sum is sparse in some—possibly distinct—basis, we can write down an EW-style set of coupled algebraic equations that use messenger fields to iteratively produce the Wiener filter solution. This requires us to add additional messenger fields. In the case where the covariance is broken into two pieces, $\mathbf{N}=\mathbf{N}_{0}+\mathbf{N}_{1}$, we can define two messenger fields ($t_{0},t_{1}$ with covariances $\mathbf{T_{0},T_{1}}$) and an auxiliary (data-like) field $d_{0}$. Then the following set of equations iteratively yields the solution to the Wiener filter equation, in analogy to the EW method:

[TABLE]

Appendix A describes the Bayesian hierarchical model that is the origin of these equations. Each additional piece of the noise covariance will describe two additional ($d_{0}$- and $t$-like) fields and equations that we can chain together to get from the data to the signal estimate. Choosing $\mathbf{T}_{0}=\tau_{0}\mathbf{I}$ and $\mathbf{T}_{1}=\tau_{1}\mathbf{I}$ makes the matrix inversions in equations (7)–(10) sparse and trivial in the proper basis, and allows iterative solutions to a much broader class of problems. The convergence proceeds fastest when these $\tau_{0},\tau_{1}$ parameters are as large as possible, and match the minimum eigenvalues of $\mathbf{N}_{0}$ and $\mathbf{N}_{1}$ respectively. (Any larger and the $\bar{\mathbf{N}}$ matrices are not positive definite.) In some cases these equations can be combined to speed the convergence (appendix A).

In our implementations, we use equations (7)–(10) directly and initialize $s=t_{0}=d_{0}=0$, but this is not critical to converge to the solution. For the problems of interest in this paper, we let $\mathbf{N}_{1}$ represent a real-space component of the noise on the full sky and $\mathbf{N}_{0}$ represent a spherical-harmonic-space component. The signal covariance is isotropic and sparse in harmonic space for all the cases we consider. We compute equation (7) in pixel space and equations (8)–(10) in harmonic space, where the matrices are diagonal and trivially invertible. Because the Wiener-filtered solution is ultimately computed in harmonic space in this implementation, it has a strict band limit.

The main computational cost is in the spherical harmonic transforms. Two transforms are required per iteration in the test problems that follow, one to take $t_{1}$ to harmonic space and one to take $d_{0}$ to real space. This is the same number of transforms required by the EW method for similar problems on the sky.

When we implement a cooling scheme, we set $\tau_{0}\rightarrow\lambda_{0}\tau_{0}$ and $\tau_{1}\rightarrow\lambda_{1}\tau_{1}$. The cooling parameters $\lambda_{0},\lambda_{1}$ start large so that the $\mathbf{T}$ matrices dominate the other covariances (the $\bar{\mathbf{N}}$ matrices are left fixed), and are gradually lowered to unity to achieve the Wiener filter solution.

We want to monitor the quality of the solution as it is converging to help us set a schedule for the cooling parameters. An obvious candidate for this job is $\chi^{2}(s)$, but unfortunately for this case, it is difficult to compute the second term in

[TABLE]

because we do not have a quick way to invert the sum of the noise matrices, as they are sparse in different bases. Conjugate gradient descent would work to compute $\chi^{2}$, but is not practical as a way to monitor the progress of our iterative solution, which purposefully avoids gradient descent techniques.

To help, we invert equation (2) to introduce a reconstructed data field, which is a function of our Wiener-filter solution and is easy to compute:

[TABLE]

In practice, we treat the harmonic-space signal and noise covariances with a band limit, so not all parts of the data are recoverable, and we do not expect to reconstruct the data exactly, particularly on small scales. Even so, the difference between the reconstructed data and the actual data,

[TABLE]

is a useful and practical metric that allows us to set the cooling schedule.

3 Results

We consider three test problems to probe the new method for Wiener filtering. This first is trivial, containing uniform white noise, and can be solved directly as well as by the iterative methods. The second has non-uniform, uncorrelated noise, which both the EW method and the new method can handle. The third addresses an application with a dense, composite noise matrix that is suitable for our extension but cannot be handled by the original EW method.

3.1 Tests with homogeneous, uncorrelated noise

In the trivial case of homogeneous noise, we can directly solve the Wiener filter (equation 3) in harmonic space, so it makes a sensible starting point to verify our Wiener filter solution. We constructed a first test problem on the sphere so that the data on large scales was dominated by an isotropic signal (with a red power spectrum), while on small scales it was dominated by uniform white noise, added in real space with covariance $\mathbf{N}=(\sigma^{2}/\Omega_{\rm pix})\mathbf{I}$. (Equivalently, the noise has power spectrum $N_{l}=\sigma^{2}$.) Then we solved for the Wiener-filtered map three ways: (1) directly in harmonic space, (2) using the EW iterative method, and (3) using our multiple messenger field method. We used a HEALPix111http://healpix.sourceforge.net pixelization at $N_{\rm side}=64$ resolution. All harmonic space computations were limited to $l\leq 3N_{\rm side}=192$.

Both iterative methods converge to the direct, harmonic space solution. For the EW method, in the limit $\tau_{\rm EW}\rightarrow\sigma^{2}$, the messenger field represents all the uniform noise, and the method reduces to the direct solution and trivially converges in a single step. For this problem the $\chi^{2}$ for the EW method is slightly better than for the harmonic space solution ($\Delta\chi^{2}/\chi^{2}\sim\mbox{few}\times 10^{-6}$), perhaps because the EW solution partially handles the noise in pixel space where the noise is generated, and the harmonic space solution includes only noise power up to the band limit.

The implementation of the multiple-messenger method converges to the same solution, but the uniform noise case is awkward for it because we need to make a decision about how to split the noise power between the $\mathbf{N}_{0}$ and $\mathbf{N}_{1}$ parts of the covariance. (The code fails if one of the covariances is left completely empty.) Putting the bulk of the white noise into $\mathbf{N}_{0}$ gives the fastest convergence: if a ninety-nine percent of the noise power goes into $\mathbf{N}_{0}$ the $\chi^{2}$ value converges to a part in $10^{3}$ after four iterations, but if only half the power is assigned to $\mathbf{N}_{0}$, this level of convergence takes 75 iterations without any cooling scheme. This method converges to a $\chi^{2}$ better than the harmonic solution, but not as good as the EW solution (and again, $\chi^{2}$ values for all solutions have only slight differences, $\Delta\chi^{2}/\chi^{2}\sim\mbox{few}\times 10^{-6}$).

3.2 Tests with inhomogeneous, uncorrelated noise

Our second test is the case that the EW solution handles most straightforwardly: inhomogeneous but uncorrelated noise. We modify the above test data set to include a mask (where the noise covariance for pixels is infinite, and the inverse covariance is zero) and to boost the noise in pixels near the mask. Both the EW method and the multiple messenger field method converge to the same solution. Empirically the convergence of the new method’s $\chi^{2}$ is exponential (as it is for the EW method). The rate of convergence depends again on the decomposition of the noise matrix and was about about 100 times slower than EW in some cases, although we did not try to optimize for speed or implement any cooling for this comparison. The noisy regions at the mask edge were the slowest to converge.

3.3 Application with dense noise

Our third test case is one that our new multiple messenger field method can handle, but the EW method cannot. We construct a noise covariance that is dense both in real and harmonic space. In this way it mimics ground-based CMB experiments that observe the statistically isotropic microwave background, but also noise that is spatially anisotropic and contains correlations from the atmosphere and scan pattern. In recent measurements from, for example, ACTPol (Louis et al., 2016), atmospheric noise dominates the temperature signal at large scales, while detector noise dominates the signal at small angular scales. At the edges of the map, where less integration time is spent, the noise level is also higher. Our test problem tries to incorporate all these features (and a mask) but treats them at reduced resolution (compared to ACTPol) to shorten runtimes for study. We use $N_{\rm side}=256$ resolution and harmonic band limit $l\leq 768$.

We build our mock data from an isotropic signal, a Gaussian random field ($s_{\rm CMB}$) generated from a $\Lambda$CDM power spectrum for CMB temperature and smoothed with a 0.3 deg full-width-half-maximum beam, plus two noise terms:

[TABLE]

The $n_{0}$ portion of the noise is a Gaussian random field, and has a covariance that is diagonal in harmonic space, but not isotropic:

[TABLE]

We show the shape of the power spectrum for the correlated noise in Figure 1, but half of the power is put into the $m=0$ modes alone, which introduces horizontal noise stripes. We also give this harmonic component of the noise half of the white noise for the cleanest portions of the map. Figure 2 shows a realization of this noise component along with the other fields.

The $n_{1}$ portion of the noise is also a Gaussian random field, but has a covariance that is diagonal is pixel space. The mask roughly corresponds to the observable sky from the Atacama, excluding the portions outside declination range $-68^{\circ}<\delta<22^{\circ}$. It also excludes the Galactic plane within $|b|<20^{\circ}$. The mask is enforced by giving masked pixels infinite variance, or in practice setting the inverse variance $\bar{N}_{1,pp}^{-1}=0$ for masked pixel index $p$. Otherwise, the noise is spatially uncorrelated but the pixel variance rises near the edges of the mask. In the least noisy portions of the sky, the signal power exceeds the noise for $139<l<588$.

For the specific example we consider, we start the cooling parameters at $\lambda\sim 5\times 10^{5}$, and iterate until the quantity $(\Delta d^{\dagger}\Delta d)$ changes by a fraction smaller than $10^{-2}$ per iteration. At that time, we reset $\lambda\rightarrow\lambda^{0.7}$, and continue to iterate. The procedure gradually lowers $\lambda$ to unity.

In Figure 2 we show the Wiener filter solution, and in Figure 3, we show how the cooling, the reconstruction’s distance to the data, and the reconstruction error each converge as we iterate. The Wiener filter solution has power only on intermediate scales, cleanly approaches zero in the masked regions, and shows little to no residual horizontal striping. With this specific cooling prescription, the $\lambda$ parameters approach unity after about 50 iterations. (We have not tried hard to optimize the cooling.) At roughly that time, the reconstructed data makes its closest approach to the actual data. It moves away upon further iterations, but the reconstruction error $\epsilon=s-s_{\rm true}$, inaccessible for real data but accessible for our test problem, continues to improve. Implementation of the cooling schedule was vital to obtain a good solution. Convergence without cooling is very much slower (constant $\lambda=1$ in Figure 3), and stripes persisted even in the masked regions after a few thousand iterations.

4 Discussion

The test problems demonstrate that our extension to the EW method can compute Wiener filtered solutions, and treat cases with complicated noise properties. We can straightforwardly extend such an approach to CMB polarization. Ground-based $B$-mode observatories will provide detailed measurements of the covariance structure of the atmosphere in polarization at microwave frequencies, allowing construction of the proper noise covariance models, although this may require significant work. Other extensions can apply to CMB lensing or large scale structure measurements that have a combination of local and correlated noise effects.

Because of the close connection of the Wiener filter both to quadratic power spectrum estimation (Oh et al., 1999) and Gibbs sampling (Wandelt et al., 2004; Eriksen et al., 2004; Jasche & Lavaux, 2015), multiple-messenger methods may allow improvements to power spectrum estimation for ground-based CMB observations. We will continue to explore these avenues at higher resolution and in more realistic scenarios applied to, for example, the Simons Observatory and CMB-S4.

The usefulness of this approach ultimately depends on our ability to construct noise covariance matrices as a sum $\mathbf{N}=\sum_{i}\mathbf{N}_{i}$ so that the inversions of $\bar{\mathbf{N}}_{i}=\mathbf{N}_{i}-\mathbf{T}_{i}$ are efficient. Already in this category are all covariances based on pixel-space and harmonic-space masks and observation hit-count maps. Parallel noise stripes, as considered here, could be placed in any orientation via straightforward rotations of the map. Small numbers of specific modes, represented as

[TABLE]

for a small matrix $\mathbf{\Lambda}$, may also be treated efficiently. Inversion via the Woodbury formula,

[TABLE]

is efficient because the term in parenthesis is a small matrix (with the dimension of $\mathbf{\Lambda}$), and the other matrix multiplications and inversions of $\mathbf{T}$ are trivial.

Finally, deliberate modifications of the covariance matrices away from their realistic values can achieve other desirable outcomes. As is standard practice, intentionally letting the covariance for specific (perhaps untrustworthy) modes go to infinity will project out those modes entirely. This strategy can downweight large-scale ground pickup or sidelobe features. For another example, we could set the signal covariance to be artificially small and diagonal, say $\mathbf{S}=\alpha\mathbf{I}$ for sufficiently small $\alpha$, to allow computation of approximately inverse noise weighted maps in cases where the noise covariance is not easily inverted. By casting it as a Wiener filtering problem,

[TABLE]

these multiple-messenger methods can find an approximate solution, given a suitable decomposition of $\mathbf{N}$.

Acknowledgments

KMH thanks Aditya Rotti for valuable discussions during the course of this work. KMH also thanks Benjamin Wandelt for providing useful comments on a draft of this paper. We acknowledge support from the NASA ATP program under grant NNX17AF87G. Some of the results in this paper have been derived using the HEALPix package (Górski et al., 2005).

Appendix A Iterative Wiener filtering for a composite noise covariance

The Wiener filter maximizes the Gaussian probability of the signal given the data, signal covariance, and noise covariance. The logarithmic probability in this case, up to a normalization constant, is

[TABLE]

We obtain the Wiener filter solution by finding the signal vector which maximizes the conditional probability for fixed data, $P(s|d)$. We proceed by completing the square on the right hand side, which yields

[TABLE]

where $\mu_{s}=(\mathbf{S}^{-1}+\mathbf{N}^{-1})^{-1}\mathbf{N}^{-1}d$. Thus the mean (maximum probability) point of the distribution, $s=\mu_{s}$, is the Wiener solution in equation (3).

We can arrive at the Elsner & Wandelt (2013, hereafter EW) solution by first noting that the augmented Gaussian probability distribution

[TABLE]

is equivalent to the above equation (19) for $\mathbf{\bar{N}=N-T}$ after the messenger field $t$ is marginalized out. In other words,

[TABLE]

The distribution is Gaussian, and so marginalizing one variable does not change the maximum probability positions for the other variables. Thus the value of $s$ that maximizes the probability is the same in both distributions $P_{1}$ and $P$.

Completing the square for $t$ and $s$ respectively yields the EW solution (our equations (5)–(6)). By iterating, we repeatedly find the maximum probability point along slices of constant $s$ and $t$, and step our way to the maximum probability for the joint distribution $P_{1}$, which in the end gives the same value for the signal $s$ as the Wiener solution from $P$.

Furthermore, Jasche & Lavaux (2015) showed that in a Bayesian hierarchical framework, this scheme trivially lends itself to Gibbs sampling. We sample the posterior probability for the signal and signal covariance $P(s,\mathbf{S}|d)$ by iteratively sampling the conditional distributions

[TABLE]

The first two conditional distributions are Gaussian and arise from completing the square, and the third is an inverse Gamma (or inverse Wishart) distribution (Larson et al., 2007).

The original EW work considered the total covariance as the sum of two parts, for the signal and the noise. The equations they derive are always valid, but we have noted that the method is really suited to the case when the noise and signal covariance are each sparse in convenient bases. Here we generalize to show that if one of the covariances is still dense, but decomposes into pieces that are sparse, we can continue to break it up until we get pieces that are easily invertible.

We proceed by adding additional fields. We treat the case with $\mathbf{N}=\mathbf{N}_{0}+\mathbf{N}_{1}$, but the idea generalizes easily if more parts are needed to represent the covariance. A hierarchical forward model in this case (Fig. 4) builds up the noise bit-by-bit. The Gaussian probability corresponding to this diagram is

[TABLE]

where $\bar{\mathbf{N}}_{0}=\mathbf{N}_{0}-\mathbf{T}_{0}$ and $\bar{\mathbf{N}}_{1}=\mathbf{N}_{1}-\mathbf{T}_{1}$.

In the same way as before, we complete the square for the signal, messenger, and auxiliary fields in (A). This provides in turn the Gaussian probability distributions for each variable conditioned on the others. The means of these distributions yield our iterative solution to the Wiener filter, equations (7)–(10). The full distributions can be used for Gibbs sampling, in analogy to equation (23).

As with EW, smart choices for the covariances of the messenger fields are proportional to the identity matrix so that they are sparse in any basis,

[TABLE]

Thus the messenger field $t_{0}$ represents $s$ plus some uniform noise; auxiliary field $d_{0}$ includes the remainder of the $\mathbf{N}_{0}$ noise; messenger field $t_{1}$ includes some more uniform noise; and $d$ includes the remainder of the $\mathbf{N}_{1}$ noise.

Because we compute both equations (8)–(9) in the sparse basis for $\mathbf{N}_{0}$, they can be combined as

[TABLE]

which solves for $t_{0}$ in terms of $t_{1}$ and $s$. Note also that equations (9)–(10) resemble the EW solution for the Wiener filter, and so

[TABLE]

useful in the specific case where the signal covariance is sparse in the same basis as the $\mathbf{N}_{0}$ portion of the noise. These shortcuts may speed up the convergence in some cases.