Issue 
A&A
Volume 550, February 2013



Article Number  A15  
Number of page(s)  10  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201220332  
Published online  21 January 2013 
Lowℓ CMB analysis and inpainting
^{1}
Laboratoire AIM, UMR CEACNRSParis 7, Irfu, Service d’Astrophysique, CEA
Saclay,
91191
GifsurYvette Cedex,
France
email: jstarck@cea.fr
^{2}
GREYC CNRSENSICAENUniversité de Caen,
6 Bd du Maréchal
Juin, 14050
Caen Cedex,
France
^{3}
Laboratoire d’Astrophysique, École Polytechnique Fédérale de
Lausanne (EPFL), Observatoire de Sauverny, 1290
Versoix,
Switzerland
Received:
4
September
2012
Accepted:
26
October
2012
Reconstructing the cosmic microwave background (CMB) in the Galactic plane is extremely difficult due to the dominant foreground emissions such as dust, freefree or synchrotron. For cosmological studies, the standard approach consists in masking this area where the reconstruction is insufficient. This leads to difficulties for the statistical analysis of the CMB map, especially at very large scales (to study for instance the low quadrupole, integrated Sachs Wolfe effect, axis of evil, etc.). We investigate how well some inpainting techniques can recover the lowℓ spherical harmonic coefficients. We introduce three new inpainting techniques based on three different kinds of priors: sparsity, energy, and isotropy, which we compare. We show that sparsity and energy priors can lead to extremely highquality reconstruction, within 1% of the cosmic variance for a mask with a sky coverage larger than 80%.
Key words: methods: data analysis / techniques: image processing / cosmic background radiation
© ESO, 2013
1. Introduction
Because component separation methods do not provide good estimates of the cosmic microwave background (CMB) in the Galactic plane and at locations of point sources, the standard approach for a CMB map analysis is to consider that the data are not reliable in these areas, and to mask them. This leads to an incomplete coverage of the sky that has to be handled properly for the subsequent analysis. This is especially true for analysis methods that operate in the spherical harmonic domain where localization is lost and fullsky coverage is assumed. For power spectrum estimations, methods such as MASTER (Hivon et al. 2002) solve a linear illposed inverse problem that allows to deconvolve the observed power spectrum of the masked map from the mask effect.
For a nonGaussianity analysis, many approaches to deal with this missing data problem have been proposed. Methods to solve this problems are called gapfilling or inpainting methods.
The simplest approach is to set the pixel values in the masked area to zero. This is sometimes claimed to be a good approach by assuming it does not add any information. However, this is not correct, as setting the masked area to zero actually adds (or removes) information, which again differs from the true CMB values. A consequence of this gapfilling technique is the creation of many artificial high amplitude coefficients at the border of the mask in a wavelet analysis or a leakage between different multipoles in a spherical harmonic analysis. This effect can be reduced with an apodized mask. A slightly more sophisticated method consists of replacing each missing pixel by the average of its neighbors and iterating until the gaps are filled. This technique is called diffuse inpainting and has been used in PlanckLFI data preprocessing Zacchei et al. (2011). It is acceptable for treatment of point sources, but is far from being a reasonable solution for the Galactic plane inpainting in a CMB nonGaussianity analysis.
In Abrial et al. (2007; 2008), the problem was considered as an illposed inverse problem, z = Mx, where x is the unknown CMB, M is the masking operator, and z the masked CMB map. A sparsity prior in the spherical harmonic domain was used to regularize the problem. This sparsitybased inpainting approach has been successfully used for two different CMB studies, the CMB weaklensing on Planck simulated data (Perotto et al. 2010; Plaszczynski et al. 2012), and the analysis of the integrated SachsWolfe effect (ISW) on WMAP data (Dupé et al. 2011). In both cases, the authors showed from MonteCarlo simulations that the statistics derived from the inpainted maps can be trusted at a high confidence level, and that sparsitybased inpainting can indeed provide an easy and effective solution to the large Galactic mask problem.
It was also shown that sparse inpainting is useful for weaklensing data (Pires et al. 2009), Fermi data (Schmitt et al. 2010), and asteroseismic data (Sato et al. 2010). The sparse inpainting success has motivated the community to investigate more inpainting techniques, and other approaches have recently been proposed.
Bucher & Louis (2012) and Kim et al. (2012) sought a solution that complies with the CMB properties, i.e. to be a homogeneous Gaussian random field with a specific power spectrum. Good results were derived, but the approach presents the drawback that we need to assume a given cosmology, which may not be appropriate for nonGaussianity studies. For large scale CMB anomaly studies, Feeney et al. (2011) and BenDavid et al. (2012), building upon the work of de OliveiraCosta & Tegmark (2006), proposed to use generalized leastsquares (GLS), which coincide with the maximumlikelihood estimator (MLE) under appropriate assumptions such as Gaussianity. This method also requires an input power spectrum. In addition, when the covariance in the GLS is singular, the authors proposed to regularize it by adding a weak perturbation term that must also to ensure positivedefiniteness. With this regularization, the estimator is no longer a GLS (nor a MLE).
Because the missing data problem is illposed, prior knowledge is needed to reduce the space of candidate solutions. Methods that handle this problem in the literature assume priors, either explicitly or implicitly. If we exclude the zeroinpainting and the diffuse inpainting methods, which are of little interest for the Galactic plane inpainting, the two other priors are:

Gaussianity: the CMB is assumed to be a homogeneous andisotropic Gaussian random field, therefore it is adequate to usesuch a prior. In practice, methods using this prior require thetheoretical power spectrum, which has either to be estimatedusing a method such as TOUSI (Paykariet al. 2012), or a cosmology has tobe assumed, which is an even stronger assumption than theGaussianity prior. We should also keep in mind that the goal ofnonGaussianity studies is to check that the observed CMB doesnot deviate from Gaussianity. Therefore, we should be carefulwith this prior.

Sparsity: the sparsity prior on a signal consists of assuming that its coefficients in a given representation domain, when sorted in decreasing order of magnitude, exhibit a fast decay rate, typically a polynomial decay with an exponent that depends on the regularity of the signal. Spherical harmonic coefficients of the CMB show a decay in O(ℓ^{2}) up to ℓ around 900 and then O(ℓ^{3}) for larger multipoles. Thus, the sparsity prior is advocated, and this explains its success for CMBrelated inverse problems such as inpainting or component separation (Bobin et al. 2012).
We here revisit the Gaussianity and sparsity priors, and introduce an additional one, the CMB isotropy, to recover the spherical harmonic coefficients at low ℓ (<10) from masked data. We describe novel and fast algorithms to solve the optimization problems corresponding to each prior. These algorithms originate from the field of nonsmooth optimization theory and are efficiently appliedto largescale data. We then show that some of these inpainting algorithms are very efficient to recover the spherical harmonic coefficients for ℓ < 10 when using the sparsity or energy priors. We also study the reconstruction quality as a function of the sky coverage, and we show that a very good reconstruction quality, within 1% of the cosmic variance, can be reached for a mask with a sky coverage better than 80%.
2. CMB inpainting
2.1. Problem statement
We observe z = Mx, where x is a realvalued centered and squareintegrable field on the unit sphere S^{2}, and M is a linear bounded masking operator. The goal is to recover x from z.
The field x can be expanded as where the complexvalued functions Y_{ℓm} are the cocalled spherical harmonics, ℓ is the multipole moment, and m is the phase ranging from − ℓ to ℓ. The a_{ℓ,m} are the spherical harmonic coefficients of x. In the following we will denote S the spherical harmonic transform operator and S^{∗} its adjoint (hence its inverse since spherical harmonics form an orthobasis of L_{2}(S^{2})).
If x is a widesense stationary (i.e., homogeneous) random field, the spherical harmonic coefficients are uncorrelated, Moreover, if the process is isotropic, then where C_{ℓ} is the angular power spectrum, which depends solely on ℓ.
In the remainder of the paper, we consider a finitedimensional setting, where the sphere is discretized. x can be rearranged in a column vector in R^{n}, and similarly for z ∈ R^{m}, with m < n, and a ∈ C^{p}. Therefore, M can be seen as a matrix taking values in { 0,1 }, i.e., M ∈ M_{m × n}({ 0,1 }). The goal of inpainting is to recover x from z.
2.2. a general inpainting framework
The recovery of x from z when m < n amounts to solving an underdetermined system of linear equations. Traditional mathematical reasoning, in fact the fundamental theorem of linear algebra, tells us not to attempt this: there are more unknowns than equations. However, if we have prior information about the underlying field, there are some rigorous results showing that such an inversion might be possible (Starck et al. 2010).
In general, we can write the inpainting problem as the following optimization program (1)where R is a proper lowerbounded penalty function reflecting some prior on x, and is closed constraint set expressing the fidelity term. Typically, in the noiseless case, . This is what we focus on in the remainder of the paper. We consider three types of priors, each corresponding to a specific choice of R.
3. Sparsity prior
Sparsitybased inpainting has been proposed for the CMB in (Abrial et al. 2007, 2008) and for weaklensing mass map reconstruction in Pires et al. (2009, 2010). In Perotto et al. (2010), Dupé et al. (2011) and Rassat et al. (2012), it was shown that sparsitybased inpainting does not destroy CMB weaklensing, ISW signals or some largescale anomalies in the CMB, and is therefore an elegant way to handle the masking problem. The masking effect can be thought of as a loss of sparsity in the spherical harmonic domain because the information required to define the map has been spread across the spherical harmonic basis (leakage effect).
Optimization problem
If the spherical harmonic coefficients a of x (i.e. a = Sx) are assumed to be sparse, then a wellknown penalty to promote sparsity is the l_{q} (pseudo or quasi)norm, with q ∈ [0,1]. Therefore, Eq. (1)becomes (2)where , and where for z ∈C. For q = 0, the l_{0} pseudonorm counts the number of nonzero entries of its argument. The inpainted map is finally reconstructed as .
Solving Eq. (2)when q = 0 is known to be NPhard. This is additionaly complicated by the nonsmooth constraint term. Iterative Hard thresholding (IHT), described in Appendix A, attempts to solve this problem. It is also wellknown that the l_{1} norm is the tightest convex relaxation (in the ℓ_{2} ball) of the l_{0} penalty (Starck et al. 2010). This suggests solving Eq. (2)with q = 1. In this case, the problem is wellposed: it has at least a minimizer, and all minimizers are global. Furthermore, although it is completely nonsmooth, it can be solved efficiently with a provably converging algorithm that belongs to the family of proximal splitting schemes (Combettes & Pesquet 2011; Starck et al. 2010). This can be done with the DouglasRachford (DR) algorithm described in Appendix B.
Comparison between IHT and DR
Fig. 1
Relative MSE per ℓ in percent for the two sparsitybased inpainting algorithms, IHT (dashed red line) and DR (continuous black line). 

Open with DEXTER 
We compared the IHT and DR sparsitybased inpainting algorithms in 100 MonteCarlo simulations using a mask with sky coverage Fsky = 77%. In all our experiments, we used 150 iterations for both iterative schemes, β = 1 and α_{n} ≡ 1 (∀n) in the DR scheme (see Appendix B). For each inpainted map i ∈ {1,··· ,100}, we computed the relative mean squarederror (MSE) and the its version in percent per ℓFigure 1 depicts the relative MSE in percent per ℓ for the two sparsitybased inpainting algorithms IHT (l_{0}) and DR (l_{1}). We see that at very low ℓ, l_{1}sparsity inpainting as provided by the DR algorithm yields better results. We performed the test with other masks and arrived at similar conclusions. Because we focus on lowℓ, only the l_{1} inpainting as solved by the DR algorithm is considered below.
4. Energy prior
Optimization problem
If we know a priori the powerspectrum (C_{ℓ,m})_{ℓ,m} (not the angular one, which implicitly assumes isotropy) of the Gaussian field x, then using a maximum a posteriori (MAP) argument, the inpainting problem amounts to minimizing a weighted ℓ_{2}norm subject to a linear constraint (3)where for a complexvalued vector a, , i.e. a weighted ℓ_{2} norm. By strong convexity, Problem Eq. (3)is wellposed and has exactly one minimizer, therefore justifying equality instead of inclusion in the argument of the minimum in Eq. (3).
Fig. 2
Set B_{ℓ}. 

Open with DEXTER 
Fig. 3
Normalized critical thresholds at the significance level 0.05 as a function of ℓ. 

Open with DEXTER 
Since the objective is differentiable with a Lipschitzcontinuous gradient, and the projector Proj_{{ x:z = Mx }} on the linear equality set is known in closedform, one can use a projected gradientscheme to solve Eq. (3). However, it turns that the estimate of the descent stepsize can be rather crude, which may hinder the convergence speed. One can think of using an inexact linesearch, but this will complicate the algorithm unnecessarily. This is why we propose a new algorithm in Appendix C, based on the DouglasRachford splitting method, which is easy to implement and efficient.
In all our experiments in the remainder of the paper, we used 150 iterations for this algorithm with β = 50 and α_{n} ≡ 1 (∀n) (see Appendix C). This approach requires the powerspectrum C_{ℓ,m} as an input parameter. In practice, an estimate of the power spectrum from the data using MASTER was used, and we set . The latter is reminiscent of an isotropy assumption, which is not necessarily true.
5. Isotropy prior
5.1. Structural constraints on the power spectrum
Strictly speaking, we would define the set of isotropic homogeneous random processes on the (discrete) sphere as the closed set where C_{ℓ,m} = E[(Sx)_{ℓ,m}^{2}] and C_{ℓ} is the angular power spectrum, which depends solely on ℓ. Given a realization x of a random field, we can then naively state the orthogonal projection of x onto by solving This formulation of the projection constraint set is not straightforward to deal with for at least the following reasons: (i) the isotropy constraint set involves the unknown true C_{ℓ} (through the expectation operator), which necessitates resorting to stochastic programming; (ii) the constraint set is not convex.
5.2. Projection with a deterministic constraint
An alternative to the constraint set would be to replace C_{ℓ,m} and C_{ℓ} by their empirical sample estimates, i.e., C_{ℓ,m} by a_{ℓ,m}^{2} and C_{ℓ} by . However, this hard constraint might be too strict in practice and we propose to make it softer by taking into account the variability inherent to the sample estimates, as we explain now.
In a nutshell, the goal is to formulate a constraint set, where for each ℓ, the entries a_{ℓ,m} of the spherical harmonic coefficient vector a deviate the least possible in magnitude (up to a certain accuracy) from some prespecified vector μ; typically we take or its empirical estimate . Put formally, this reads is a compact set, although not convex.
Fig. 4
Six masks with Fsky ∈ { 0.98,0.93,0.87,0.77,0.67,0.57 }. The masked area is in blue. 

Open with DEXTER 
Fig. 5
Top: inputsmoothed simulated CMB map (ℓ_{max} = 10) and the same map, but masked (Fsky = 77%) and not smoothed (i.e., inputsimulated data). Middle: inpainting of the top right image up to ℓ_{max} = 10, using a sparsity prior (left) and an energy prior (right). Bottom: inpainting of the top right image up to ℓ_{max} = 10, using an isotropy prior (left) and a simple Fsky correction (right). 

Open with DEXTER 
Fig. 6
ℓ_{1} sparsitybased inpainting of a simulated CMB map with different masks. Top: Fsky = 98% and 93%, middle: 87% and 77%, and bottom: 67% and 57%. 

Open with DEXTER 
We now turn to the projection on . We begin by noting that this set is separable, i.e., ,
where B_{ℓ} is the band depicted Fig. 2. It turns out that for fixed μ, B_{ℓ} is socalled proxregular since its associated orthogonal projector is singlevalued with a closedform. Indeed, the projector onto is (4)where (5)
Choice of the constraint radius:
to devise a meaningful choice (from a statistical perspective) of the constraint radius ε, we first need to derive the null distribution^{1} of a_{ℓ,m} − μ_{ℓ},∀(ℓ,m), with the particular case when . a_{ℓ,m} being the spherical harmonic coefficients of a zeromean stationary and isotropic Gaussian process whose angular power spectrum is C_{ℓ}, it is easy to show that where L = 2ℓ + 1. By the Fisher asymptotic formula, we then obtain where means convergence in distribution. Thus, ignoring the obvious dependency between a_{ℓ,m} and μ_{ℓ}, we obtain The distribution of this difference can be derived properly using the delta method in law, but computing the covariance matrix of the augmented vector (a_{ℓ,m})_{m}, remains a problem. The quality of the above asymptotic approximation increases as ℓ increases.
The upper and lower critical thresholds at the doublesided significance level 0 ≤ α ≤ 1 are where Φ is the standard normal cumulative distribution function.
We depict in Fig. 3 the upper and lower critical thresholds normalized by at the classical significance level 0.05 as a function of ℓ. One can observe that the thresholds are not symmetric. This entails that different values of ε should be used in Eq. (5). Furthermore, as expected, the two thresholds decrease in magnitude as ℓ increases. They attain a plateau for ℓ that is sufficiently high (typically ≥ 100). It is easy to see that the two limit values are and .
The isotropyconstrained inpainting algorithm is given is Appendix D.
6. Experiments
6.1. a_{ℓ,m} reconstruction
In this section, we compare the different inpainting methods described above: the DRsparsity prior inpainting, the DRenergy prior inpainting, and the isotropy prior inpainting. We also test the case of noinpainting by applying just an Fsky correction to the a_{ℓ,m} spherical harmonic coefficients of the map (i.e., data a_{ℓ,m} are corrected with a correction term equal to ). We ran 100 CMB MonteCarlo simulations with a resolution corresponding to nside equal to 32, which is precise enough for the large scales we are considering. We considered six different masks with sky coverage Fsky of { 0.98,0.93,0.87,0.77,0.67,0.57 }.
Figure 4 displays these masks. Point source masks were not considered here, since they should not affect a_{ℓ,m} estimation at very low multipole. Each of these simulated CMB maps was masked with the six masks. Figure 5 top right shows one CMB realization masked with the 77% Fsky mask. Figure 5 middle and bottom show the inpainting of the Fig. 5 top right image with the four methods (i.e., using sparsity, energy, isotropy priors, and Fsky correction). Figure 6 depicts the results given by the ℓ_{1} sparsitybased inpainting method for the six different masks. We can clearly see that the quality of the reconstruction degrades with decreasing sky coverage.
Fig. 7
Relative MSE (in percent) per ℓ versus sky coverage for different multipoles. See text for details. 

Open with DEXTER 
The a_{ℓ,m} coefficients of the six hundred maps (100 realizations masked with the six different masks) were then estimated using the three inpainting methods and the Fsky correction methods. The approaches were compared in terms of the relative MSE and relative MSE per ℓ; see end of Sect. 3 for their definition.
Figure 7 shows the relative MSE for ℓ = 2 to 5 versus the sky coverage Fsky. The three horizontal lines correspond to 1%, 5%, and 10% of the cosmic variance (CV). Inpainting based on the sparsity and energy priors give relatively close results, which are better than the one assuming the isotropy prior or the Fsky correction. It is interesting to notice that a very high quality reconstruction (with 1% of the CV) can be obtained with both the sparsity and energybased inpainting methods up to ℓ = 4 for a mask with an Fsky exceeding 80%. If WMAP data do not allow us to study nonGaussianity with such a small mask, Planck component separation will certainly be able to achieve the required quality to enable make possible the useusing of such a small Galactic mask. With a 77%coverage Galactic mask, the error increases by a factor 5 at ℓ = 4 !
Figure 8 shows the same errors, but they are now plotted versus the multipoles for the six different masks. The ℓ_{1} sparsitybased inpainting method seems to be slightly better than the one based on the energy prior, especially when the sky coverage decreases.
Fig. 8
Relative MSE in percent per sky coverage versus ℓ for the six masks. 

Open with DEXTER 
6.2. Anomalies in CMB maps
One of the main motivations of recovering fullsky maps is to be able to study statistical properties such as Gaussianity and statistical isotropy on large scales. Statistical isotropy is violated if there exists a preferred axis in the map. Mirror parity, i.e., parity with respect to reflections through a plane: , where is the normal vector to the plane, is an example of a statistic where a preferred axis can be sought (BenDavid et al. 2012; Land & Magueijo 2005).
With fullsky data, one can estimate the Smap for a given multipole in spherical harmonics by considering (BenDavid et al. 2012): (6)where corresponds to the value of the a_{ℓm} coefficients when the map is rotated to have as the zaxis. Positive (negative) values of correspond to even (odd) mirror parities in the direction. The same statistic can also be considered summed over all low multipoles one wishes to consider(e.g., focusing only on low multipoles as in BenDavid et al. 2012): (7)It is convenient to redefine the parity estimator as , so that ⟨S⟩ = 0.
The most even and odd mirrorparity directions for a given map can be considered by estimating (BenDavid et al. 2012) (8)(9)where μ(S) and σ(S) are the mean and standard deviation of the S map. The quantities S_{+} and S_{−} each correspond to different axes and and are the quantities which we consider when discussing mirror parity in CMB maps.
To use inpainted maps to study mirror parity, it is crucial that the inpainting method constitutes a biasfree reconstruction method, i.e., that it does not destroy existing mirror parities nor that it creates previously inexistent mirror parity anomalies.
In order to test this, we considered three sets of 1000 simulated CMB maps using a WMAP7 bestfit cosmology with nside = 32, one set for each inpainting prior. For each simulation we calculated the mirror parity estimators S_{±}. We isolated the anomalous maps, defined as those whose S_{±} value is higher than twice the standard deviation given by the simulations. This yielded 34 maps (3.4%) for the even mirror parity and 40 maps (4%) for the odd mirror parity in the fullsky simulations.
Fig. 9
Percentage difference between anomalous (i.e., high) odd mirror parity scores (S_{−}) before (i.e., fullsky) and after inpainting for three different inpainting priors. 

Open with DEXTER 
Since recent work showed a potential oddmirror anomaly in WMAP data (BenDavid et al. 2012), we focused on testing for potential biases in oddmirror anomalies. Figure 9 shows the percentage difference between the true S_{−} value and the one estimated from the inpainted map. For this statistical test, the result does not depend on the inpainting method because the three different priors give similar biases and error bars. For positivemirror anomalies, we find similar results (not shown in figure). Similarly to the previous experiment, we can see that a very good estimation of the parity statistic can be achieved for masks with Fsky exceeding 0.8.
7. Software and reproducible research
To support reproducible research, the developed IDL code will be released with the next version of ISAP (Interactive Sparse astronomical data Analysis Packages) via the web site: http://www.cosmostat.org
All experiments were performed using the default options

Sparsity:
Alm = cmb_lowl_alm_inpainting(CMBMap, Mask, /sparsity)

Energy:
Pd = mrs_powspec( CMBMap, Mask)
P = mrs_deconv_powspec( pd, Mask)
Alm = cmb_lowl_alm_inpainting(CMBMap, /Energy, Prea=P)

Isotropy:
Pd = mrs_powspec( CMBMap, Mask)
P = mrs_deconv_powspec( Pd, Mask)
Alm = cmb_lowl_alm_inpainting(CMBMap, /Isotropy, Prea=P)
8. Conclusion
We have investigated three priors to regularize the CMB inpainting problem: Gaussianity, sparsity, and isotropy. To solve the corresponding minimization problems, we also proposed fast novel algorithms, based for instance on proximal splitting methods for convex nonsmooth optimization. We found that both Gaussianity and ℓ_{1} sparsity priors lead to very good results. The isotropy prior does not provide comparably good results. The sparsity prior seems to lead to slightly better results for ℓ > 2, and more robust when the sky coverage decreases. Furthermore, unlike the energyprior based inpainting, the sparsitybased one does not require a power spectrum as an input, which is a significant advantage.
Then we evaluated the reconstruction quality as a function of the sky coverage and saw that a highquality reconstruction, within 1% of the cosmic variance, can be reached for a mask with an Fsky exceeding 80%. We also studied mirrorparity anomalies and found that mask with an Fsky exceeding 80% also permitted nearly biasfree reconstructions of the mirror parity scores. Such a mask seems unrealistic for a WMAP data analysis, but it could and should be a target for Planck component separation methods. To know if this 80% mask is realistic for Planck, we will need more realistic simulations including the ISW effect, and residual foregrounds at amplitudes similar to those achieved by the actual best Planck component separation methods. Another interesting study will consist of comparing the sparse inpainting to other methods such as the maximum likelihood or the Wiener filtering (BenDavid et al. 2012; Feeney et al. 2011).
Appendix A: Algorithm for the l _{0} problem
Solving Eq. (2)when q = 0 is known to be NPhard. This is further complicated by the presence of the nonsmooth constraint term. Iterative Hard thresholding (IHT) attempts to solve this problem through the following scheme (10)where the nonlinear operator is a termbyterm hard thresholding, i.e. , where if  a_{i} > λ and 0 otherwise. The threshold parameter λ_{n} decreases with the iteration number and is supposed to play a role similar to the cooling parameter in simulated annealing, i.e. hopefully it allows to avoid stationary points and local minima.
Appendix B: Algorithm for the l _{1} problem
It is now wellknown that l_{1} norm is the tightest convex relaxation (in the ℓ_{2} ball) of the l_{0} penalty (Starck et al. 2010). This suggests solving Eq. (2)with q = 1. In this case, the problem is wellposed: it has at least a minimizer, and all minimizers are global. Furthermore, although it is completely nonsmooth, it can be solved efficiently with a provably convergent algorithm that belongs to the family of proximal splitting schemes (Combettes & Pesquet 2011; Starck et al. 2010).
In particular, we propose to use the DouglasRachford (DR) splitting scheme. Let β > 0, be a sequence in ]0, 2[ such that ∑ _{t ∈N}α_{n}(2 − α_{n}) = +∞. The DR recursion, applied to Eq. (2)with q = 1, reads (11)prox_{β∥·∥1}is the proximity operator of the l_{1}norm, which can be easily shown to be softthresholding (12)where and . Proj_{{ a:z = MS∗a }} is the orthogonal projector on the corresponding set, and it consists of taking the inverse spherical harmonic transform of its argument, setting its pixel values to the observed ones at the corresponding locations, and taking the forward spherical harmonic transform. It can be shown that the sequence converges to a global minimizer of Eq. (2)for q = 1.
Appendix C: Algorithm for the l _{2} problem,
This scheme is based again on DouglasRachford splitting applied to solve problem Eq. (3). Its steps are (13)where β and α_{n} are defined as above, and the proximity operator of the squared weighted ℓ_{2}norm is (14)where ⊗ stands for the entrywise multiplication between two vectors, and the division is also to be understood entrywise. The projector Proj_{{x:z = Mx}} has a simple closedform and consists in setting pixel values to the observed ones at the corresponding
locations, and keeping the others intact. It can be shown that the sequence converges to the unique global minimizer of Eq. (3).
Appendix D: Algorithm for isotropy inpainting
The isotropyconstrained inpainting problem can be cast as the nonconvex feasibility problem (15)The feasible set is nonempty since the constraint set is nonempty and closed; the CMB is in it under the isotropy hypothesis. To solve Eq. (15), we propose to use the von Neumann method of alternating projections onto the two constraint sets, whose recursion can be written (16)Using closedness and proxregularity of the constraints, and by arguments from Lewis & Malick (2008), we can conclude that our nonconvex alternating projections algorithm for inpainting converges locally to a point of the intersection which is nonempty.
Acknowledgments
The authors would like to thank Hiranya Peiris and Stephen Feeney for useful discussions. This work was supported by the French National Agency for Research (ANR 08EMER00901), the European Research Council grant SparseAstro (ERC228261), and the Swiss National Science Foundation (SNSF).
References
 Abrial, P., Moudden, Y., Starck, J., et al. 2007, J. Fourier Anal. Appl., 13, 729 [CrossRef] [Google Scholar]
 Abrial, P., Moudden, Y., Starck, J., et al. 2008, Stat. Meth., 5, 289 [NASA ADS] [CrossRef] [Google Scholar]
 BenDavid, A., Kovetz, E. D., & Itzhaki, N. 2012, ApJ, 748, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Bobin, J., Starck, J.L., Sureau, F., & Basak, S. 2012, A&A, in press DOI: 10.1051/00046361/201219781 [Google Scholar]
 Bucher, M., & Louis, T. 2012, MNRAS, 3271 [Google Scholar]
 Combettes, P. L., & Pesquet, J.C. 2011, FixedPoint Algorithms for Inverse Problems in Science and Engineering (SpringerVerlag), 185 [Google Scholar]
 de OliveiraCosta, A., & Tegmark, M. 2006, Phys. Rev. D, 74, 023005 [NASA ADS] [CrossRef] [Google Scholar]
 Dupé, F.X., Rassat, A., Starck, J.L., & Fadili, M. J. 2011, A&A, 534, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feeney, S. M., Peiris, H. V., & Pontzen, A. 2011, Phys. Rev. D, 84, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., Naselsky, P., & Mandolesi, N. 2012, ApJ, 750, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Land, K., & Magueijo, J. 2005, Phys. Rev. Lett., 72, 101302 [Google Scholar]
 Lewis, A., & Malick, J. 2008, Math. Oper. Res., 33, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Paykari, P., Starck, J.L., & Fadili, M. J. 2012, A&A, 541, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perotto, L., Bobin, J., Plaszczynski, S., Starck, J., & Lavabre, A. 2010, A&A, 519, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pires, S., Starck, J., Amara, A., et al. 2009, MNRAS, 395, 1265 [NASA ADS] [CrossRef] [Google Scholar]
 Pires, S., Starck, J., & Refregier, A. 2010, IEEE Sign. Proc. Mag., 27, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Plaszczynski, S., Lavabre, A., Perotto, L., & Starck, J.L. 2012, A&A, 544, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rassat, A., Starck, J.L., & Dupe, F. X. 2012, A&A, submitted [Google Scholar]
 Sato, K. H., Garcia, R. A., Pires, S., et al. 2010, Astron. Nachr., submitted [arXiv:1003.5178] [Google Scholar]
 Schmitt, J., Starck, J. L., Casandjian, J. M., Fadili, J., & Grenier, I. 2010, A&A, 517, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J.L., Murtagh, F., & Fadili, M. 2010, Sparse Image and Signal Processing (Cambridge University Press) [Google Scholar]
 Zacchei, A., Maino, D., Baccigalupi, C., et al. 2011, A&A, 536, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1
Relative MSE per ℓ in percent for the two sparsitybased inpainting algorithms, IHT (dashed red line) and DR (continuous black line). 

Open with DEXTER  
In the text 
Fig. 2
Set B_{ℓ}. 

Open with DEXTER  
In the text 
Fig. 3
Normalized critical thresholds at the significance level 0.05 as a function of ℓ. 

Open with DEXTER  
In the text 
Fig. 4
Six masks with Fsky ∈ { 0.98,0.93,0.87,0.77,0.67,0.57 }. The masked area is in blue. 

Open with DEXTER  
In the text 
Fig. 5
Top: inputsmoothed simulated CMB map (ℓ_{max} = 10) and the same map, but masked (Fsky = 77%) and not smoothed (i.e., inputsimulated data). Middle: inpainting of the top right image up to ℓ_{max} = 10, using a sparsity prior (left) and an energy prior (right). Bottom: inpainting of the top right image up to ℓ_{max} = 10, using an isotropy prior (left) and a simple Fsky correction (right). 

Open with DEXTER  
In the text 
Fig. 6
ℓ_{1} sparsitybased inpainting of a simulated CMB map with different masks. Top: Fsky = 98% and 93%, middle: 87% and 77%, and bottom: 67% and 57%. 

Open with DEXTER  
In the text 
Fig. 7
Relative MSE (in percent) per ℓ versus sky coverage for different multipoles. See text for details. 

Open with DEXTER  
In the text 
Fig. 8
Relative MSE in percent per sky coverage versus ℓ for the six masks. 

Open with DEXTER  
In the text 
Fig. 9
Percentage difference between anomalous (i.e., high) odd mirror parity scores (S_{−}) before (i.e., fullsky) and after inpainting for three different inpainting priors. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.