Free Access
Issue
A&A
Volume 566, June 2014
Article Number A100
Number of page(s) 12
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201322706
Published online 23 June 2014

© ESO, 2014

1. Introduction

Missions such as Planck and WMAP provide high-resolution full-sky maps of microwave emissions. The primary goal of these missions is to measure the fluctuations in the cosmic microwave background (CMB), the relic radiation of the big bang, and hence provide key information about the birth and evolution of our Universe. These missions also provide full-sky information on the Galactic and extragalactic emissions in a frequency range that has not been probed before. One of the main scientific challenges in the CMB data analysis consists of accurately separating the CMB, the various Galactic emissions (synchrotron, free-free, dust and Galactic compact-source emissions, to name a few) as well as the contribution from extragalactic sources. In particular, bright compact-source emissions cover a significant fraction of the sky (Planck Collaboration VII 2011; Wright et al. 2009), even at high Galactic latitudes, where the CMB is less affected by other diffuse foregrounds. These emissions therefore need to be carefully dealt with for CMB or Galactic data analysis.

From a source-separation perspective, these emissions are difficult to model since the compact sources display both spectral and temporal variability. As reported in Wright et al. (2009), at least one-third of the extragalactic WMAP sources display temporal variability with high confidence, with more than a 2:1 range in fluxes. On the other hand, extrapolating the fluxes from catalogues obtained at lower frequency than that probed by WMAP or Planck is erroneous for radio sources (Planck Collaboration XIV 2011): flat-spectrum emissions from the core of extragalactic sources are dominating for Planck or WMAP channels, whereas steep-spectrum lobe emissions dominate at lower frequencies. Finally, their high spectral variability makes it very difficult to estimate their flux with generic source-separation techniques that implicitly rely on the factorization of spatial and spectral information (Tegmark & de Oliveira-Costa 1998; Vielva et al. 2003; López-Caniego et al. 2006). This means that each CMB data set needs to be processed independently to accurately estimate the compact-source fluxes in this data set.

The most distinctive information for separating compact and diffuse emissions is based on morphology, in particular for point sources where the shape of the emission is the point spread function (PSF) of the instrument at the wavelength considered. Consequently, specific pipelines were developed based on this information in the CMB data analysis; the compact sources are first detected and flux-estimated and then are either masked or their contribution subtracted before component separation is performed (Planck Collaboration XII 2014).

Numerous approaches based on the morphology of the sources have been proposed to address the initial problem of source detection in the CMB data: morphological image-processing such as Sextractor (Bertin & Arnouts 1996), matched filters (Tegmark & de Oliveira-Costa 1998; Barreiro et al. 2003; López-Caniego et al. 2006), wavelet-based techniques (Cayón et al. 2000; Vielva et al. 2001; Martínez-González et al. 2003; Barreiro et al. 2003; López-Caniego et al. 2006), and Bayesian detection (Hobson & McLachlan 2003; Savage & Oliver 2007; Carvalho et al. 2009; Guglielmetti et al. 2009; Argüeso et al. 2011). Multichannel techniques, which take into account either an estimate of the power spectrum of diffuse emissions (Herranz & Sanz 2008) or the alleged spectral signatures for the point sources (Lanz et al. 2011; Vio et al. 2013) have also been devised. Other multichannel techniques detect the point sources by cancelling the CMB contribution in the data by a weighted sum of the channels (Wright et al. 2009; Ramos et al. 2011; Scodeller et al. 2012). We refer to Paykari & Starck (2012) for a brief review of the techniques. Some compact-source catalogues have been published for the WMAP (for instance, Wright et al. 2009; Chen & Wright 2009; Ramos et al. 2011; Scodeller et al. 2012; Bennett et al. 2013) or Planck (Planck Collaboration VII 2011; Planck Collaboration XXVIII 2014) missions.

When bright sources are detected, their contribution can be removed from the data. Several strategies have been proposed.

  • Masking: this is a standard approach where a point-source mask isapplied to the data map before it is analysed. This solution is veryefficient for CMB power spectrum estimation, and the price to payis only a small loss of sensitivity at low multipoles (Noltaet al. 2009; Larsonet al. 2011). However, for CMBnon-Gaussianity analysis the mask may create spurious featuresin the non-Gaussianity estimators. In addition, maskingcomplicates the Galactic diffuse emission analysis

  • Masking and inpainting: to avoid artefacts caused by the mask on large scales, the masked data can be interpolated using inpainting techniques. Several methods have been proposed, such as diffuse inpainting (Zacchei et al. 2011), CMB-constrained realization inpainting (Bucher & Louis 2012; Kim et al. 2012), or sparse inpainting (Abrial et al. 2007, 2008; Starck et al. 2013b). It was shown that sparse inpainting, based on the sparsity of the CMB in spherical harmonics, does not significantly impact non-Gaussianity measures such as the skewness and kurtosis of the CMB data (Abrial et al. 2007) or the integrated Sachs-Wolfe signal (Dupé et al. 2011) and the weak-lensing signal (Perotto et al. 2010). These inpainting approaches are, however, expected to be inaccurate in estimating the CMB in the masked regions, which would also be detrimental to Galactic studies.

  • Fitting: a more progressive approach is to first model the background due to diffuse emissions, then model the shape of the compact sources by modelling the instrument PSF, and finally estimate the flux by aperture photometry or PSF fitting. This is typically performed when the flux and the position are not jointly estimated in the detection algorithm, as performed in the Bayesian detection approaches (where the background can be modelled through a covariance matrix or described by a low-order polynomial). This method was adopted by the Planck and WMAP teams to estimate the point-source flux (Planck Collaboration VII 2011; Planck Collaboration XXVIII 2014; Wright et al. 2009; Bennett et al. 2013), where the background was modelled as a low-order polynomial (baseline or linear background). The few background parameters and the point-source flux were then estimated locally by χ2 minimization.

Scodeller & Hansen (2012) showed that removing the point-source fluxes instead of masking them provides similar results on the final CMB power spectrum estimation with no biases on higher order statistics. A cleaned CMB map is much easier to analyse than a masked CMB map, since there is no need to handle the point-source mask anymore. Another advantage is that it provides a CMB estimate at the point-source positions in the mask (or a point-source-free channel estimate if applied on a channel instead of on the reconstructed CMB map), useful in particular for Galactic studies, while inpainting techniques only try to fill the gaps without destroying the statistical properties of the map.

In this paper, we propose a new approach for point-source removal that is based on a morphological source-separation method, assuming the sources have already been detected. We use a more flexible (and potentially more complex) model for the background to capture its fluctuations more accurately at the scale of the sources. In this approach we assume that the diffuse emissions are sparse in the spherical harmonic domain, while the point sources are sparse in the direct domain, and resolved compact sources are sparse in the wavelet domain. In Sect. 2, we describe this method along with the proposed algorithm, which is adapted from recent convex-optimization techniques that solve the corresponding inverse problem. Results on compact-source removal in full-sky WMAP realistic simulations are presented in Sect. 3, where our approach is compared with the standard flux-fitting or low-order polynomial background-fitting as performed by WMAP and Planck consortia (Planck Collaboration VII 2011; Planck Collaboration XXVIII 2014; Wright et al. 2009; Bennett et al. 2013). Sparse point-source removal is then applied on WMAP data in Sect. 4, code information is given in Sect. 5, and conclusions and perspectives are drawn in Sect. 6.

2. Sparse point-source removal

2.1. Modelling of the sky

We modelled microwave-channel data over the full sky y ∈ RNe (i.e. the channel map is composed of Ne real-valued pixels) composed of three components { x1,x2,x3 } ∈ RNe: the point sources, the extended compact sources, and the diffuse emissions. In this decomposition, the diffuse background is composed of the CMB itself and the synchrotron, free-free, dust diffuse emissions as well as all the sources with flux below the detection cut. The forward model can be cast as (1)where n ∈ RNe is an additional noise, assumed to be a real-valued Gaussian random field (but not necessarily stationary). In the associated inverse problem, we need to estimate three unknowns from one equation, which is not possible without additional constraints/information. We now discuss these assumptions and constraints on each of the components.

  • Point-source catalogues are available for WMAP and Planckdata sets. Using these catalogues, x1 can be modelled as a set of Dirac functions { δp } p convolved by the instrumental beam b: where Np is the number of point sources that need to be removed and fp is the flux of the pth point source. The source fluxes are expected and enforced to be positive.

  • In case of extended compact sources, the morphological information is weaker. Several extended sources were detected and flagged in the Planck data set, mostly located close to the Galactic centre, while in the WMAP data such sources were not considered (detection was performed outside a Galactic mask). In this work, we assumed that these emissions can be modelled through a sparse decomposition in the undecimated isotropic spherical wavelet dictionary (Starck et al. 2006): , where Ns is the number of wavelet scales, θjk corresponds to a wavelet atom at scale j and position k, and we have overall Nw = NeNs wavelet coefficients at all scales. Assuming sparsity of this Galactic component, only a few wavelet coefficients wjk at selected locations k close to the Galactic centre are assumed significantly away from zero.

  • The background emission x3 is much more complex to model accurately. It is usually represented as a local low-order polynomial background (either baseline, or first order) because of its local smoothness, and in this work it is characterized by its sparsity in the spherical harmonic domain. This assumption of sparsity in CMB applications is motivated in Starck et al. (2013a). x3 is also assumed band-limited because of the beam effect.

These components can only be separated if the components are sparse enough in mutually incoherent dictionaries, such as a dirac basis and a spherical harmonic basis (Donoho & Huo 2001; Donoho & Elad 2003). Said differently, compact-source emissions are expected to decrease the sparsity level of the data measured on spherical harmonics, which will be the key ingredient driving the separation process in our approach.

2.2. Associated inverse problem

Using all these hypotheses and constraints, we consider the following inverse problem: (2)where

  • f ∈ RNp is a vector containing the flux of the Np point sources, is the positive orthant of RNp to enforce the positivity constraint, and B ∈ RNe × Np is a matrix operating on fluxes and implementing the local projection of the beam at the location of the point sources;

  • is a vector containing the Nw wavelet coefficients at all scales, W represents the matrix implementing the spherical wavelet transform (containing each wavelet atom as a column), M corresponds to a binary mask equal to one where compact extended sources are expected. We promote sparsity of x2 by making use of the 1 norm of w as a penalization term;

  • a ∈ CNl contains the Nl complex spherical harmonic coefficients of x3 and is the space of band-limited signals on the sphere (up to a chosen multipole max depending on the resolution of the channel). S ∈ CNe × Nl is the matrix describing the orthogonal spherical harmonic transform: x3 = Sa. Sparsity of x3 is enforced through the use of a 1 norm penalization term on a, which is computed as the sum of the modulus of the complex multipole as in the sparse inpainting described in Starck et al. (2013a)

  • Σ is the (non-stationary) noise covariance matrix and denotes the square of the 2 norm weighted by Σ-1.

The reconstructed point-source-free map can be obtained by . The problem described in Eq. (2) is a convex problem and is related to a constrained morphological component analysis (Starck et al. 2004, 2005) or a basis pursuit denoising problem (Chen & Donoho 1998) with a deconvolution step. The principal advantage of this constrained formulation lies in having only a few hyperparameters to set, which are easily interpretable: the expected noise level ϵ; the balance between 1 norm in the wavelet and spherical harmonic domains. More explicitely, the higher the ratio (γ/β), the sparser the solution for x3 is in the spherical harmonic domain, and conversely, the lower the ratio, the sparser the estimate of x2 in the wavelet domain. Setting them equal would lead to penalize both non-sparse solutions for x2 in the wavelet domain and x3 in the spherical harmonic domain, as performed in morphological component analysis (Starck et al. 2004, 2005). Note that once the ratio is set, changing their modulus does not change the inverse problem.

The main difficulty here lies in controlling the interplay between sparsity and data fidelity constraints: how can we efficiently estimate a sparse solution (more precisely with minimal 1 norm), knowing that many combinations of x1, x2 and x3 can satisfy the data fidelity constraint?

In the next section, we propose an algorithm to solve this problem.

2.3. Proposed algorithm

Inverse problems as described by Eq. (2) can be solved using primal-dual approaches (e.g. Chambolle & Pock 2011; Briceño Arias & Combettes 2011; Becker et al. 2011; Combettes & Pesquet 2012). The proposed algorithm was derived from Chambolle & Pock (2011), it requires only one application of the costly spherical harmonic and wavelet transforms and one application of their adjoint per iteration and does not require sub-iterations. The algorithm is set as follows:

morphological component analysis with a primal-dual approach (SPSR)

  • 1-

    Choose (a0,f0,w0,t0) ∈ CNl × RNp × RNw × RNe and , , . Choose also the hyperparameters γ,β and τ,σ s.t. (assuming normalized | | S | | 2 = | | MW | | 2 = | | B | | 2 = 1).

  • 2-

    Iterate (n ≥ 0): (3)

where is the projection onto the positive orthant (i.e. sets to 0 negative fluxes) and is the projection onto the set of considered band-limited signals (i.e. sets to 0 all multipoles greater than the chosen max); is the standard soft-thresholding operator applied component-wise: (4)where | xi | is the complex modulus of xi for complex vectors. From Chambolle & Pock (2011), the sequence ((fn,an,wn),tn) converges to a saddle point of the primal-dual problem with a restricted dual gap decreasing as (first-order method). We initialised the algorithm with null images.

The sparse-source-removal (SPSR) algorithm has four parameters to set: β, γ, τ, and σ. The last two are the primal and dual steps, respectively, and drive the convergence speed of the algorithm. The ratio of the first two parameters is the hyperparameter that controls the balance in the sparsity of the decomposition of the extended compact sources and the diffuse background as discussed in the previous section. The choice of their modulus also affects the convergence speed of the algorithm: a high modulus leads to a slowly built sparse approximation, or said differently, the algorithm would provide very sparse solutions at the beginning of the algorithm that are far from satisfying the data fidelity constraint, however; a low modulus, in contrast, leads to building approximations that satisfy the data fidelity constraint early in the iterations, but are not as sparse. In practice, several values need to be experimentally tested to obtain an algorithm with reasonable convergence speed and to derive the required number of iterations.

To assess the relative performance of the proposed approach, flux estimates obtained with SPSR were compared with the low-order polynomial fitting approach (Bennett et al. 2013; Planck Collaboration VII 2011; Planck Collaboration XXVIII 2014). Fluxes were estimated in a 3.5σ region as previously recommended in Wright et al. (2009) using a Levenberg-Marquardt algorithm with local χ2 minimization (using the C++ library ALGLIB (www.alglib.net), Sergey Bochkanov and Vladimir Bystritsky). The local background was either modelled with a baseline or a first-order polynomial (FIT-C and FIT-L in the following) and a positivity constraint was applied on the fluxes. Clustered point sources with overlapping fitting regions were jointly fitted for. However, in practice, such situations were scarce in the WMAP simulations reported.

3. Results on synthetic WMAP simulations

3.1. Planck Sky Model simulations at WMAP frequencies

The Planck Sky Model (PSM) software (Delabrouille et al. 2013) was employed to simulate WMAP-like data. Each one of the five WMAP frequency channels included the following:

  • A diffuse component, comprised of a Gaussian CMB, generatedfrom a six-parameter Λ-CDM model (with default values from WMAP 7-year data combined with BAO and the Hubble constant measurements, and with added C lensing) and synchrotron, free-free, thermal-dust and spinning-dust emissions, simulated with default PSM parameters.

  • A compact component, constituted by radio, infrared, and strong ultra-compact HII region emissions, as well as a far-infrared background (model jgn2005 in the PSM). More specifically, radio sources are essentially derived from observations conducted between 0.85 GHz and 4.85 GHz. Each one of these ~2 million simulated sources possessed its own spectrum, composed of power laws with spectral indices varying across pre-defined bands of electromagnetic frequencies. As reported in Delabrouille et al. (2013), this catalogue is expected to faithfully reproduce the clustering properties of the radio sources and the observed/modelled source number counts at 5 and 20 GHz. In the PSM, the infrared sources were derived from the IRAS point-source catalogue and faint-source catalogue, with assumed modified blackbody emissions to extrapolate the flux at lower frequencies. Simulated sources were also added near the Galactic plane and in IRAS gaps to obtain the same mean surface density as a function of flux (down to 80 mJy) as in the regions well covered by IRAS. Emissions from ultra-compact HII regions were also derived from IRAS, with flux extrapolated according to a modified blackbody fit with added low-frequency flux with a free-free spectral index when a radio-counterpart is found. Note that compact emissions from the Galaxy were also simulated as part of diffuse emissions. Finally, the simulated far-infrared background is composed of realistic distribution of clustered point sources with assumed flux-dependent spectral indices and reproduced Planck observations reasonably well.

  • Noise modelled as a non-stationary Gaussian random field, with variance in each pixel derived from the WMAP 9-year hit maps; the channel bandpasses were modelled as diracs located at the centre frequency, and the beams were assumed Gaussian with a full-width at half maximum given by the beam size provided by the WMAP consortium.

Overall, the PSM simulations reflect the complexity in estimating the strong point-source fluxes (which is essential to test the robustness of the proposed approach), with complex diffuse background, unresolved sources with clustering properties below the detection limit, strong compact emissions in the Galactic plane, and non-stationary noise. A more comprehensive description of the PSM is provided in Delabrouille et al. (2013). A patch extracted from these WMAP 9-year simulated data is represented in Fig. 1 at each of the five WMAP frequency channels to illustrate the complexity of the source separation, because the resolution, noise, and background change in each of the channels. At low-frequency channels, the background fluctuates at the scale of the point sources, and noise contribution to the flux estimates are lower than in the high-frequency channels, where noise becomes important at the scale of the source, as reported in Wright et al. (2009).

thumbnail Fig. 1

Patch of the simulated sky centred on two detected point sources at the WMAP wavelengths.

Open with DEXTER

thumbnail Fig. 2

Evolution with the number of iterations of the normalized 2 norm of the residuals (left) and normalized 1 norm of the diffuse and compact component for channel Ka and for various values for the modulus of β = γ. Reaching the upper limit of the data-fidelity constraint leads to a normalized residual 2 norm of 1. The 1 norm was normalized to the 1 norm of the input data in direct space (sum of the absolute values of the pixels in the frequency map). M corresponds to the highest modulus of the spherical harmonic coefficients in the input image.

Open with DEXTER

3.2. Flux recovery in noise-limited catalogue

3.2.1. Catalogue generation

First, a catalogue of point sources was independently generated from the simulations for each WMAP channel by retaining the sources with flux above 5σ of the estimated radiometer noise (Wright et al. 2009). The number of sources retained for each channel and the flux cut are displayed in Table 1. This scenario reflects an optimistic situation where only noise would be a limiting factor for source detection (and not the background). Consequently, it allows us to investigate how well point sources can be estimated even when their flux is similar to the background levels. Note that the cut in flux increases with frequency (because of the increasing level of the noise) and our sources are mainly radio sources (i.e. with flux decreasing with increasing frequency). Consequently, fewer sources are detected at higher frequency channels.

Table 1

Characteristics of the point-source catalogue generated from simulations.

thumbnail Fig. 3

Flux recovery with the proposed approach as a function of the number of iterations. After 13 350 (1350) iterations, the flux estimates have converged for channel K (W). Note that the fluxes are displayed in logarithmic scales to visualize the convergence speed for both low- and high-flux sources.

Open with DEXTER

Table 2

Statistics on the recovered flux for each channel for the three approaches and for various flux bands (in mJy).

3.2.2. SPSR working conditions

The working conditions of SPSR were chosen as follows: the noise hyperparameter ϵ corresponds in the following to the 95th percentile of a χ2(Ne) distribution according to the whitened noise statistics. For all channels the hyperparameter ratio was set to β/γ = 1 to favour both sparsity of x2 in the wavelet domain and x3 in the spherical harmonic domain, as discussed before.

A range of values for the modulus of β = γ was first tested and turned out to be crucial for the convergence speed of the algorithm. This is illustrated in Fig. 2, where SPSR was run with a high number of iterations for channel Ka. Setting the modulus too high (or too low) leads to many iterations necessary to fulfill the data-fidelity constraint (or minimizing the 1 penalty for components that satisfy the data-fidelity term). For instance, for values of β = 0.001M (or β = 0.005M) the data-fidelity constraint is satisfied in the first one thousand iterations, but the sparsity-penalty term decreases very slowly and has not converged after 50 000 iterations. Conversely, β = 0.1M (or β = 0.05M) leads to very sparse approximations that do not satisfy the data-fidelity constraint even after 50 000 iterations. We chose the value β = 0.01M for all channels, which led to the best convergence properties among the tested values.

thumbnail Fig. 4

Statistics on error flux using the various approaches for the five channels for five different flux bands (flux <500 mJy, 500 mJy flux <1 Jy, 1 mJy flux <2 Jy, 2 Jy flux <5 Jy, flux 5 Jy). Quartiles and extreme values are plotted.

Open with DEXTER

thumbnail Fig. 5

Point-source residuals with noise for the various approaches in a projected region of channel Q, illustrating the lower bias obtained for this channel with the proposed approach. These maps were obtained by subtracting the diffuse component from the point-source-subtracted map. In the first two maps, three negative regions can be visually detected because of the overfitting of the point-source flux in the FIT-C and FIT-L approaches, which is different from the SPSR residual map. This does not display biases that are significantly higher than the noise level.

Open with DEXTER

Finally, the number of iterations needs to be set. In addition to the modulus of β = γ, the convergence speed of the algorithm depends on the conditioning of the beam matrix B. The number of iterations therefore varies from channel to channel, requiring more iterations for lower resolution. The flux estimates at several stages of the proposed algorithm are illustrated in Fig. 3 for channels K and W, which correspond to the lowest and highest resolution and to the highest and lowest noise, respectively. In practice, 13 350 iterations were chosen for channel K and 9750 iterations for all other channels as a compromise between processing time and level of convergence. We also ensured that the statistics measured were marginally changing with respect to the relative difference between methods when reaching this number of iterations. Generally, the large number of iterations in both cases is caused by the absence of regularity in the inverse problem of Eq. (2), which leads to first-order algorithms such as SPSR. To give a pratical example of the computational cost of SPSR, 9750 iterations for the channel Ka represent ~80 h of processing in a cluster with 20 Intel Xeon 2 GHz processors, essentially spent in the spherical harmonic and wavelet transforms.

thumbnail Fig. 6

Projected regions with large discrepancy between the proposed approach and the local low-order polynomial background modelling. Images correspond to the simulated data, the diffuse component with noise, and the point-source-subtracted data with the three approaches. SPSR in this case allows a better subtraction of the point-source contribution because the local variations in the background are not captured by the other two approaches.

Open with DEXTER

3.2.3. Flux recovery and point-source subtraction

We compared point-source flux recovery in different flux bands for the FIT-C, FIT-L, and SPSR approaches in the simulated channels. The bias, root-mean square errors (RMSE), and mean absolute deviation (MAD) (less sensitive to outliers) are displayed in Table 2. For all methods, the RMSE or MAD varies from channels K to W as a consequence of increasing noise contribution and decreasing background fluctuations at the scale of the point sources, as previously observed in Wright et al. (2009). When the noise contribution to flux uncertainty is predominant (in channel W), FIT-C and FIT-L have generally lower biases, RMSE, and MAD than SPSR. Note, however, that the improvement is very small (at most 6% difference in RMSE, bottom row of Table 2). Therefore, it seems that in the noise-dominating case, fitting a few local parameters for the diffuse background, as in FIT-C and FIT-L, gives slightly more robust flux estimates than SPSR, where many more global parameters are estimated. The situation is the opposite when background fluctuations are the foremost contribution as in the first four channels: SPSR now leads to much lower RMSE and MAD than FIT-C or FIT-L, as well as lower biases (for the first three channels). In particular, point sources with flux lower than 1 Jy are better recovered by about 12%, 26%, 32% and 8% in terms of RMSE for channels K to V, with much lower biases for the first three channels (more than 70 mJy better for channel K, for instance).

Statistics on the distribution of the error in the different flux bands are displayed in Fig. 4. In these box-plots both quartiles (three horizontal bars making the “box”) and extreme values (upper and lower horizontal bars) are represented for the various approaches and the different channels. First, the same conclusions can be drawn as in Table 1 by focusing on the interquartile range and median for the different channels: SPSR outperforms FIT-C and FIT-L, except for channel W, where results are only slightly degraded. In particular, Fig. 5 illustrates that the bias in the lowest frequency channels is higher than the noise level in the residual maps for FIT-C and FIT-L, but not for SPSR.

These box-plots also illustrate that the proposed approach leads to lower extreme errors, while FIT-L or FIT-C can both fail with large flux errors. An example of such an incorrect flux-fitting is illustrated in Fig. 6 for channel K, where FIT-C and FIT-L lead to overestimated flux in a complex background region that cannot be accurately modelled with low-order polynomials. In contrast, SPSR successfully estimates the background and therefore gives a better flux estimate.

3.3. Flux recovery with internal multichannel catalogue derived from simulations

The noise-free catalogue is useful to assess whether the more flexible background model in SPSR indeed better captures the background fluctuation than low-order local polynomials, in particular, to recover low-flux sources. However, in practice detection is performed internally for each channel, and a catalogue is built by combining these multichannel detections. The flux-recovery problem is slightly different in this case, since the catalogue is first subject to Eddington bias (excess detection in positive fluctuations of the background, which translates into higher biases in estimating the flux and impacts low-flux sources), and then the task is also to recover as accurately as possible the spectra of sources that can be lower than the detection level in some of the channels (e.g. in higher frequencies for radio sources).

3.3.1. Catalogue generation

As was previously done for the WMAP official point-source catalogue, we therefore constructed a catalogue by first applying a matched filter to each channel, taking into account second-order statistics of CMB and noise in spherical harmonic space (Bennett et al. 2013; Tegmark & de Oliveira-Costa 1998). Sources were considered detected if the value in the filtered map was higher than 5σ of the filtered background (computed locally). This resulted in 689 (442, 415, 248, and 172) sources for channel K (Ka, Q, V, and W). Compared with the results obtained in Table 1, a sub-population of the previous catalogue is essentially considered at this step: sources with sufficiently high local contrast to the diffuse background that they can be internally detected. After merging of the sources, the final simulation-derived multichannel catalogue contained 724 sources, and we processed these in the same way as desribed in the previous paragraph. However, after merging, the sources considered for a channel with this approach are no longer a sub-population of the previous catalogue: sources detected in other bands are considered as well.

Table 3

Statistics on the recovered flux for each channel using the multichannel internal catalogue for the three approaches and for various flux bands (in mJy).

thumbnail Fig. 7

Statistics on error flux computed from internal catalogue and using the various approaches for the five channels for five different flux bands (flux <500 mJy, 500 mJy flux <1 Jy, 1 mJy flux <2 Jy, 2 Jy flux <5 Jy, flux ≥ 5 Jy). Quartiles and extreme values are plotted.

Open with DEXTER

3.3.2. Flux recovery and point-source subtraction

The same working conditions were set for SPSR as in the previous simulations. The statistics obtained as in Table 2 are presented in Table 3. The values reported for medium- or high-flux point sources in the first channels are similar if not identical, because this sub-population did not change from the previous catalogue. However, significant discrepancies from Table 2 can be noted. First, the biases have increased for all methods for low-flux sources, an indication for Eddington bias. Note, however, that the bias has increased less for SPSR than for FIT-C or FIT-L in absolute values, leading to even larger differences between the methods (SPSR leads to the lowest biases for all channels, with a bias now up to more than 100 mJy lower for sources lower than 1 Jy in channel K). This is also true for the RMSE of these low-flux sources, which essentially did not change for SPSR and was slightly degraded for FIT-C or FIT-L. Note that most of the low-flux sources in the catalogue are not detected in the W band, which might explain why SPSR gave better RMSE for medium-flux sources than FIT-C or FIT-L, while this was not the case in the noise-limited catalogue. The proposed approach in this multichannel scenario therefore leads to the lowest errors, even for flux estimated in channel W.

thumbnail Fig. 8

Projected region in channel V of the WMAP 9-year data showing the largest discrepancy between SPSR and FIT-C or FIT-L estimated fluxes for two fitted sources. Images correspond to the WMAP data, the point-source-subtracted data with FIT-C, FIT-L, and SPSR approaches, and point-source-estimated data using SPSR. The differences here between the approaches are a fraction of noise and the diffuse component levels.

Open with DEXTER

The distribution of errors is reported in Fig. 7. SPSR seems as robust to large flux errors in the internal catalogue as in the noise-limited case with respect to FIT-C or FIT-L. Note the distribution of errors for channel W for low-flux sources (lower than 500 mJy), which were indeed not detected in this channel (compare with Fig. 4). Large flux errors are observed for all methods, but are less pronounced for SPSR than for FIT-C or FIT-L.

As a summary for the internal catalogue, the most striking difference compared with the noise-limited catalogue occurs in low-flux sources that are now subject to Eddington bias, whereas for high-flux the results are essentially similar. Better relative performance in this situation in terms of bias and RMSE or MAD may be attributed to two different phenomena: the proposed approach would be less affected by Eddington bias and/or lead to better estimates for sources not detected in the band.

4. WMAP 9-year processing

The proposed approach was then applied to subtract point-source emission from WMAP 9-year data. The maps were processed as follows: the differential assemblies data deconvolved with the asymmetrical part of the beam were averaged to obtain 9-year frequency band data (Bennett et al. 2013). The beam considered for each channel was obtained as the mean over all axisymetric beams provided for the differential assemblies at that frequency. A catalogue was then built by merging the two catalogues provided by the WMAP collaboration, resulting in 628 considered point sources.

The same processing (with the same working conditions) as for the simulations was then performed for both SPSR and the local low-order polynomial minimization. In the absence of ground truth, it is obviously difficult to quantitatively assess the relative performance of the algorithm as in the simulation. Only extreme cases that visually illustrate the performance of the approaches are therefore presented in Figs. 8 and 9.

Figure 8 represents a patch centred on the source with the strongest difference between SPSR and FIT-L for channel V. Because the difference between the approaches is only a fraction of the background fluctuations at larger scales than the sources, visually no difference can be observed for that channel. However, Fig. 9 illustrates an extreme case for channel K, suggesting some over-estimation of the flux in a complex background scenario with FIT-C and FIT-L, but not with SPSR, as reported in Fig. 6 in the simulation.

thumbnail Fig. 9

Projected region in WMAP 9-year data showing the large discrepancy in between SPSR and FIT-C or FIT-L estimated fluxes for two fitted sources. Images correspond to the WMAP data, the point-source-subtracted data with FIT-C, FIT-L, and SPSR approaches. The positive ring around two the fitted regions as well as a negative peak inside these regions seem to indicate that FIT-C or FIT-L are here over-fitting. The same phenomenon is absent for SPSR.

Open with DEXTER

5. Reproducible research

Table 4

List of products made available in this paper in the spirit of reproducible research, available at http://www.cosmostat.org/Products.

In the spirit of participating in reproducible research, we make all codes and resulting products that constitute the main results of this paper public. In Table 4 we list all products that are made freely available as a result of this paper, and which are available at http://www.cosmostat.org/Products.

For the five WMAP 9-year channels, SPSR was applied using the following command line in the open source package iSAP software1: >map_psfree =mrs_sparse_pointsource_removal(Map, GalMask, BeamInfo, StdMap, Niter=Niter).

For Ka to W maps, Niter was fixed to 9750, and to 13 350 for the K channel. Galmask is the WMAP point-source catalogue mask2, StdMap is the map containing the standard deviation of the noise per pixel, and BeamInfo is a structure containing the WMAP 9-year beam at each point-source position. The code WMAP9_DATA_routines.pro can be used to obtain the required maps and beams and to call the mrs_sparse_ pointsource_removal routine to derive the final products, as illustrated in the script wmap9_remove_point_sources.pro.

6. Conclusions

We proposed a new approach for detected point-source flux estimation and subtraction. Compared with the standard approach, which estimates point-source flux according to local low-order polynomial models of the background, the proposed technique is based on a global modelling of the background, which is assumed to be sparse in spherical harmonics to better capture its fluctuations. Bright point-source emissions decrease the sparsity of this background, which is the key phenomenon driving the separation process. An algorithm was adapted from recent convex optimization developments to solve the corresponding inverse problem.

We evaluated the proposed technique as well as techniques used in WMAP and Planck collaboration on realistic simulations of the WMAP microwave sky. In a noise-limited catalogue, except for channel W, where noise leads to faint differences between the estimates and a slightly poorer estimate for the proposed approach for low-flux sources (about 6%), our approach out-performs local polynomial fitting. In the internally derived catalogue, SPSR also consistently leads to

  • the lowest biases in the first three channels (up to100 mJy lower bias for sources <1 Jy in channel K)

  • the lowest RMSE in all channels and all flux bands (with values decreased by at least 5%, 25%, 28%,7%, and 3% and up to 35%, 28%, 32%, 9%, and 14% for channels K to W)

  • more robust point-source subtraction in complex background, as illustrated in specific examples

  • consequently better multichannel estimates of the point-source fluxes, which is useful to estimate more accurate spectra.

This technique was finally applied to the WMAP 9-year data deconvolved with the beam asymmetries and the resulting point-source-subtracted maps are available on-line. We focused in this paper on point sources. Even if the proposed model for extended compact sources is rather simple in SPSR, it is sufficient to derive a background accurate enough to improve on point-source flux recovery compared with FIT-C or FIT-L. However, better modelling of the extended compact sources is still required if they need be subtracted from the data. We left this task for future work.


Acknowledgments

This work was funded by the European Research Council grant SparseAstro (ERC-228261), and the Swiss National Science Foundation (SNSF). We used the Healpix3 software (Górski et al. 2005), the iSAP1 software, and WMAP data4.

References

All Tables

Table 1

Characteristics of the point-source catalogue generated from simulations.

Table 2

Statistics on the recovered flux for each channel for the three approaches and for various flux bands (in mJy).

Table 3

Statistics on the recovered flux for each channel using the multichannel internal catalogue for the three approaches and for various flux bands (in mJy).

Table 4

List of products made available in this paper in the spirit of reproducible research, available at http://www.cosmostat.org/Products.

All Figures

thumbnail Fig. 1

Patch of the simulated sky centred on two detected point sources at the WMAP wavelengths.

Open with DEXTER
In the text
thumbnail Fig. 2

Evolution with the number of iterations of the normalized 2 norm of the residuals (left) and normalized 1 norm of the diffuse and compact component for channel Ka and for various values for the modulus of β = γ. Reaching the upper limit of the data-fidelity constraint leads to a normalized residual 2 norm of 1. The 1 norm was normalized to the 1 norm of the input data in direct space (sum of the absolute values of the pixels in the frequency map). M corresponds to the highest modulus of the spherical harmonic coefficients in the input image.

Open with DEXTER
In the text
thumbnail Fig. 3

Flux recovery with the proposed approach as a function of the number of iterations. After 13 350 (1350) iterations, the flux estimates have converged for channel K (W). Note that the fluxes are displayed in logarithmic scales to visualize the convergence speed for both low- and high-flux sources.

Open with DEXTER
In the text
thumbnail Fig. 4

Statistics on error flux using the various approaches for the five channels for five different flux bands (flux <500 mJy, 500 mJy flux <1 Jy, 1 mJy flux <2 Jy, 2 Jy flux <5 Jy, flux 5 Jy). Quartiles and extreme values are plotted.

Open with DEXTER
In the text
thumbnail Fig. 5

Point-source residuals with noise for the various approaches in a projected region of channel Q, illustrating the lower bias obtained for this channel with the proposed approach. These maps were obtained by subtracting the diffuse component from the point-source-subtracted map. In the first two maps, three negative regions can be visually detected because of the overfitting of the point-source flux in the FIT-C and FIT-L approaches, which is different from the SPSR residual map. This does not display biases that are significantly higher than the noise level.

Open with DEXTER
In the text
thumbnail Fig. 6

Projected regions with large discrepancy between the proposed approach and the local low-order polynomial background modelling. Images correspond to the simulated data, the diffuse component with noise, and the point-source-subtracted data with the three approaches. SPSR in this case allows a better subtraction of the point-source contribution because the local variations in the background are not captured by the other two approaches.

Open with DEXTER
In the text
thumbnail Fig. 7

Statistics on error flux computed from internal catalogue and using the various approaches for the five channels for five different flux bands (flux <500 mJy, 500 mJy flux <1 Jy, 1 mJy flux <2 Jy, 2 Jy flux <5 Jy, flux ≥ 5 Jy). Quartiles and extreme values are plotted.

Open with DEXTER
In the text
thumbnail Fig. 8

Projected region in channel V of the WMAP 9-year data showing the largest discrepancy between SPSR and FIT-C or FIT-L estimated fluxes for two fitted sources. Images correspond to the WMAP data, the point-source-subtracted data with FIT-C, FIT-L, and SPSR approaches, and point-source-estimated data using SPSR. The differences here between the approaches are a fraction of noise and the diffuse component levels.

Open with DEXTER
In the text
thumbnail Fig. 9

Projected region in WMAP 9-year data showing the large discrepancy in between SPSR and FIT-C or FIT-L estimated fluxes for two fitted sources. Images correspond to the WMAP data, the point-source-subtracted data with FIT-C, FIT-L, and SPSR approaches. The positive ring around two the fitted regions as well as a negative peak inside these regions seem to indicate that FIT-C or FIT-L are here over-fitting. The same phenomenon is absent for SPSR.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.