Mon.Not.R.Astron.Soc.000,1–??(2010) Printed10January2011 (MNLATEXstylefilev2.2) Joint Bayesian separation and restoration of CMB from convolutional mixtures 1 K. Kayabol1,3⋆, J. L. Sanz2⋆, D. Herranz2⋆, E. E. Kuruoglu1⋆ and E. Salerno1⋆ 1 1ISTI, CNR, viaG. Moruzzi 1, 56124, Pisa, Italy 0 2IFCA, University of Cantabria, Avda. Los Castros s/n 39005, Santander, Spain 2 3Project-Team Ariana, INRIA, 2004 route des Lucioles, 06902 Sophia Antipolis, France n a J Accepted 1988December 15.Received1988December14;inoriginalform2010July16 7 ] ABSTRACT M We propose a Bayesianapproachto joint source separationand restorationfor astro- I physical diffuse sources. We constitute a prior statistical model for the source images h. by using their gradient maps. We assume a t-distribution for the gradient maps in p different directions, because it is able to fit both smooth and sparse data. A Monte - Carlo technique, called Langevin sampler, is used to estimate the source images and o all the model parameters are estimated by using deterministic techniques. r t s Keywords: Bayesiansourceseparation,astrophysicalimages,studenttdistribution, a Langevin. [ 1 v 7 1 INTRODUCTION studies (Bonaldi et al. 2007; Ricciardi et al. 2010), we as- 9 3 sume that the non-linear parameters of the mixing matrix Inferring the CMB radiation map is an important task to 1 are known with an error. Under this assumption, we re- estimatethecosmological parameters.Theforegroundradi- . construct the source maps in the pixel domain by using 1 ation contamination at related observation frequencies, the a Monte Carlo technique that has been recently developed 0 noise degradation of the instruments and the blur caused and tested on the astrophysical source separation problem 1 by the antenna apertures make this task very difficult. Un- 1 (Kayabolet al.2010).Ourmethodisanextendedversionof der Independent Component Analysis (ICA) framework, : the method in (Kayabolet al. 2010) to convolutional mix- v separation of the CMB radiation among the others has tureproblemandhasalsotheabilitytoestimatethemixing i been done by Maino et al. (2002). In (Cardoso et al. 2002; X matrix. Bedini et al. 2005; Bonaldi et al. 2007), the noise has been r taken into consideration to find the separation matrix and Studiesonseparationofconvolutionalorblurredimage a mixtures can be found in the image processing literature. the maps are obtained by using generalized Least Square Castella & Pesquet (2004) extended the contrast function (LS)solution.Wilson et al.(2008),Eriksen et al.(2008)and basedICAtechniqueinthecaseofblurring.Anthoine(2005) Kayabol et al. (2009) haveused Bayesian approach for sep- proposed to solve the same problem by adapting the exist- aration and noise removal of the maps. The point spread ing variational and statistical methods and modeling the functions of the antennas are included in Bedini & Salerno components in wavelet domain. Tonazzini & Gerace (2005) (2007) and Ricciardi et al. (2010) to estimate a paramet- use the Markov Random Field (MRF) based image prior ric mixing matrix, but they are not considered in the map in the Bayesian framework. Shwartz et al. (2008) address reconstruction process. a solution to separation of defocus blurred reflections in In this study, we focus on the problem of multi- the natural scenes by using the sparsity of the Short Time channel source separation and restoration from multi- FourierTransform (STFT)coefficientsaspriors.Inarecent channelblurredandnoisyobservationswithchannel-variant study(Tonazzini et al. 2010), multi-channelseparation and point spread functions (psf). The resolutions of the ob- deconvolution is proposed for document images. We use a served channel maps are generally different, since theaper- Bayesian formulation to include the effects of the antenna ture of the telescope beam depends on frequency. We per- apertures and solve the deblurring and map reconstruction form the source separation, the de-noising and the de- problem jointly. Since the psf’s of the antennas are known, blurring processes together. By considering the previous weeasilydefineourlikelihoodfunctionbyresortingtothem. In a Bayesian framework, we define prior densities for ⋆ E-mail: [email protected] (KK); [email protected] the source maps. Because of the blur and the noise, the (JLS); [email protected] (DH); [email protected] reconstruction problem is very badly conditioned. It means (EEK);[email protected](ES) thatwehavealreadylostsomedetailinformationontheob- 2 Kayabol et al. served image. The lost information in such a case is found pixel index. We assume that the observed images, y , k k ∈ in the high frequency contents of the images. While choos- 1,2,...,K , are some linear combinations of source im- { } ing our image prior, we consider this situation and define ages, s, l 1,2,...,L . Taking into account the effect of l ∈ { } a prior that models the distribution of the high frequency thetelescope, the observation model can bewritten as components of the image. We use the most basic high fre- L quencycomponentsoftheimage,namelyimagedifferentials. y =h a s +n (1) k k k,l l k Weobtaintheimagedifferentialsbyapplyingasimplehori- ∗ zontal and verticalgradient operator. The intensities of the Xl=1 image differentials are very sparse and have a heavy-tailed wheretheasteriskmeansconvolution,andhkisthechannel- distribution. We exploit the t-distribution as a statistical variant telescope point spread function (psf) in the k’th modelfortheimagedifferentials.Thefirstexamplesofuseof observation channel here assumed as Gaussian and circu- thet-distributionininverseimaging problemscanbefound larly symmetric. The observation model is not an instanta- in (Higdon 1994; Prudyuset al. 2001). In (Prudyuset al. neouslinearmixing,sincehk changesforeachchannel.The 2001), it is reported that the t-distribution approximates vector nk represents an iid zero-mean Gaussian noise with accuratelythewaveletcoefficientsofanimage.Inrecentpa- Σ=σk2IN covariancematrixwhereIN isanidentitymatrix. pers, it has been used for image restoration (Chantas et al. Althoughthenoiseisnothomogeneousintheastrophysical 2008) and deconvolution (Tzikas et al. 2009). maps, we assume that the noise variance is homogeneous In(Kayabol et al.2010),itisempiricallyshownthatthe within each sky patch and is also known. image differentialsoftheCMB, synchrotronanddustmaps canbemodeledbyt-distribution,whichcanbethenusedin Bayesian source separation. Since the CMB is assumed to 3 BAYESIAN FORMULATION OF beaGaussian random field,theimage differentialsofCMB ASTROPHYSICAL COMPONENT is more smooth than the other components. Its differential SEPARATION mightbemodelledasaGaussian,butinthisstudywemodel itasat-distributionbyusingthefactthatthet-distribution 3.1 Source Model approachestoaGaussian, ifitsdegreeof freedom (dof) pa- We used the self similarity based image model previously rameter goes to infinity. In computer experiments, we can proposedin(Kayabol et al.2010).Inthismodel,weassume deal with infinity byreplacing it bylarge numbers. thattheintensitiesoftheneighboringpixelsareclosedeach Using thestatistics of theimage differentials as a prior other.Toexpressapixelbyusingitsneighbors,wewritean for separation and reconstruction does not introduce new auto-regressive sourcemodelusing thefirst orderneighbors information into thedata, butemphasizes some part of the of thepixel in thedirection d: datatohelpthesolutionoftheproblem.Theimportantpart oftheBayesianimagereconstructionproblemisdetermining sl =αl,dGdsl+tl,d (2) the contribution of the prior to the solution. If it is defined where the maximum number of first order neighbors is 8 bytheuser,theexpectationsoftheusermightbeintroduced but we use only 4 neighbors, d 1,...,4 , in the main intothesolution. Itcan beuseful for natural, photographic ∈ { } vertical and horizontal directions. The matrix G is a lin- and medical images that are enhanced by the user, but in d earone-pixelshiftoperator, α istheregression coefficient the case of astrophysical images, since some of the physical l,d andtheregressionerrort isaniidt-distributedzero-mean parameters will be estimated after reconstruction, the con- l,d vector with dof parameter β and scale parameters δ . tribution of the prior must be controlled automatically. In l,d l,d To penalize thelarge regression error occurred in thesharp this study, the dof parameter of the t-distribution controls edge regions of the image, we use the t-distribution. Gen- thecontributionofthepriortothesolution,andweestimate erally in real images, except the Gaussian distributed ones, this parameter from data along the iterations. The disper- theregression errorisbettermodelledbysomeheavy-tailed sion(scale)parameterofthet-distributionisalsoestimated distribution. The t-distribution can also model the Gaus- in thealgorithm. siandistributeddata.Thereforeitisaconvenientmodelfor The organization of the paper is as follows. We intro- datawhosedistributionrangesfromCauchytoGaussian.In ducetheastrophysicalcomponentseparationprobleminthe (Kayabolet al. 2010), t-distribution has been fitted to sim- case of convolutional mixtures in Section 2. In Section 3, ulated CMB, synchrotron and dust maps and gives better we define formally the source separation problem in the resultsinthesenseof meansquareerrorwhencompared to Bayesian context, and outline the source model, the likeli- Gaussian and Cauchy densities. The multivariate probabil- hoodandtheposteriors.Thedetailsofthesourcemapsand itydensityfunctionofanimagemodelledbyat-distribution parameters estimations are given in Section 4. A number with mean µ (α )=α G s, scale δ and dof β can of simulation cases including for five different sky patches l,d l,d l,d d l l,d l,d bedefined as are given in Section 5, and finally conclusions are drawn in Section 6. Γ((N+β )/2) (s µ ,δ ,β ) = l,d T l| l,d l,d l,d Γ(β /2)(πβ δ )N/2 l,d l,d l,d 2 COMPONENT SEPARATION PROBLEM: 1+ ||sl−µl,d||2 −(N+βl,d)/(23) CONVOLUTIONAL MIXTURES CASE × β δ (cid:20) l,d l,d (cid:21) Let the kth observed pixel be denoted by y , where where Γ(.) is the Gamma function. Using a latent variable, k,i i 1,2,...,N represents the lexicographically ordered i.e. ν , the t-distribution can also be written in implicit l,d ∈ { } Joint Bayesian separation and restoration 3 form usingaGaussian andaGammadensity(Liu & Rubin p(δ y ,s ,A,Θ ) p(s Θ) l,d| 1:K 1:L −δl,d ∝ l| 1995): p(s y ,s ,A,Θ) p(y s ,A)p(s Θ) l| 1:K (1:L)−l ∝ 1:K| 1:L l| (s µ ,δ ,β )= (4) T l| l,d l,d l,d where ”–variable” expressions in the subscripts denote the removal of that variable from the variable set. δ I β β Z N(cid:18)sl|µl,d, lν,dl,dN(cid:19)G(cid:16)νl,d| 2l,d, 2l,d(cid:17)dνl,d. (5) usingTthheeMEMLemsteitmhaotdio(nLiouf&theRpuabrianm19et9e5r)siαsgl,idv,eβnl,idnaSnedctδiol,nd Weusetherepresentationin(5)toestimatetheparam- 4.3. To estimate the source images, we use a version of the eters using EM method. posterior p(sl .) augmented by auxiliary variables and find | We can write the density of s by using the image dif- the estimate by means of a Langevin sampler. The details l ferentials indifferentdirections,byassuming directionalin- are given in Section 4. dependence,as 4 p(s Θ)= (s µ (α ),δ ,β ). (6) 4 ESTIMATION OF ASTROPHYSICAL MAPS l| T l| l,d l,d l,d l,d AND PARAMETERS d=1 Y whereΘ= α ,β ,δ isthesetofallparam- Inthissection,wegivetheestimationofthemixingmatrix, 1:L,1:4 1:L,1:4 1:L,1:4 { } eter. source maps and their parameters. Weassume uniform priorsfor α andδ anduse un- l,d l,d informative Jeffrey’s prior for β ; β 1/β . l,d l,d l,d ∼ 4.1 Mixing Matrix WeassumethatthepriorofAisuniformbetween0and . 3.2 Likelihood ∞ From (11), it can be seen that the posterior density of a k,l Since the observation noise is assumed to be independent onlydependsontheGaussianlikelihood in(7).Wecanfind and identically distributed zero-mean Gaussian at each themaximum likelihood estimate of a as k,l pixel, thelikelihood is expressed as L 1 a = sTHT(y H a s )u(a ) (12) p(y s ,A) K exp W(s y ,A,σ2) (7) k,l sTl HTkHksl l k k− ki=X1,i6=l k,i i k,l 1:K| 1:L ∝ − 1:L| k k where u(a ) is the unit step function. kY=1 (cid:8) (cid:9) k,l (y H L a s) 2 W(s y ,A,σ2) = || k− k l=1 k,l l || (8) 1:L| k k 2σ2 4.2 Astrophysical Map Estimation Pk where y1:K and s1:L represent the set of all observed and We simulate the astrophysical maps from their posteriors sourceimages.ThemixingmatrixAcontainsallthemixing usinganMCMCscheme.IntheclassicalMCMCschemes,a coefficients ak,l introduced in (1). We assume uniform pri- random walk process is used to produce the proposal sam- ors for ak,l. Matrix Hk is the Toeplitz convolution matrix ples.Althoughrandom walk is simple,it affects theconver- constituted by hk introduced in (1). gence time adversely. The random walk process uses only the previous sample for producing a new proposal. Instead ofarandomwalk,weusetheLangevinstochasticequation, 3.3 Posteriors which exploits the gradient information of theenergy func- By taking intoaccount theparameters of thesource priors, tion to produce a new proposal. Since the gradient directs we write the joint posterior density of all unknownsas: the proposed samples towards the mode, the final sample setwillmostlycomefromaroundthemodeoftheposterior. p(s ,A,Θy ) p(y s ,A)p(s ,A,Θ) (9) 1:L | 1:K ∝ 1:K| 1:L 1:L The Langevin equation used in this studyis written as twhheerjoeinpt(yp1r:iKor|sd1:eLn,sAit)yiosftuhneklnikoewlinhso.oTdhaenjdoinpt(sp1:rLio,rAc,aΘn)bies skl+1=skl − 12Dg(sk1:L)+D21wl (13) factorizedasp(s α ,β ,δ )p(A)p(β ) p(δ1:L,1:4)p(α1:L1,:1L:4|).1F:Lu,1r:t4her1m:Lo,1r:e4,si1n:Lce,1:t4hesources1a:rLe,1a:4s- where g(sk1:L) = [∇slE(s1:L)]s1:L=sk1:L, ∇sl is the gradient sumed to be independent, the joint probability density of withrespecttosl andthediagonal matrixD12 containsthe thesources is also factorized as discrete time steps τ , n = 1 : N. The total energy func- l,n L tionE(s1:L)isproportionaltothenegativelogarithmofthe p(s1:L|Θ)= p(sl|Θ) (10) pthoestdeirffioursiaosn−coloegffipc(iseln|yt1is:KD,s(1:L=)−τl,2A.,HΘe)r.eF,omratthreixitDhpisixreel-, Yl=1 n,n l,n ferredtoasthediffusionmatrix.Wedetermineit bytaking Forestimating alloftheunknowns,wewritetheircon- theinverseofthediagonaloftheHessianmatrixofE(s ). ditional posteriors as 1:L Rather than the expectation of the inverse of Hessian ma- trix,weuseitsdiagonal calculated bythevalueof s at the l p(a y ,s ,A ,Θ) p(y s ,A) discrete time k as (Becker& LeCun 1989; Kayabol et al. k,l| 1:K 1:L −ak,l ∝ 1:K| 1:L p(α y ,s ,A,Θ ) p(s Θ) 2010) l,d| 1:K 1:L −αl,d ∝ l| p(β y ,s ,A,Θ ) p(s Θ)p(β ) (11) D=2[ (sk) ]−1. (14) l,d| 1:K 1:L −βl,d ∝ l| l,d hH l i 4 Kayabol et al. where (sk) = diag (s) and the operator diag . In the M (maximization) step, (19) is maximized with extractHs thlemain dia{gHonall}osfl=thskleHessian matrix. {} respecttoΘ.Tomaximizethisfunction,wealternateamong Since the random variables for the image pixel inten- thevariables αl,d, βl,d and δl,d. The solutions are found as sities are produced in parallel by using (13), the proce- sTGTs dure is faster than the random walk process adopted in α = l d l (21) l,d sTGTG s (Kayabol et al. 2009). Equation (13) produces a candidate l d d l map sample by taking into account the noise, the channel- φ (s,α ) δ = ν d l l,d (22) variantblurandthemixingmatrix.UnlikeLSsolution,this l,d h l,di N equationdoesnotcontainanymatrixinversion.Thederiva- The maximization with respect to β does not have a tion details of the equation can be found in (Kayabol et al. l,d simplesolution.Itcanbesolvedbysettingitsfirstderivative 2010). to zero: After the sample production process, the samples are applied to a Metropolis-Hastings (Hastings 1970) scheme ψ (β /2)+logβ + logν ν +1=0 (23) 1 l,d l,d l,d l,d pixel-by-pixel. The acceptance probability of any proposed − h i−h i sample is definedas min{ϕ(skl,+n1,skl,n),1}, where wdihgearmemψa1(f.u)niscttihone.first derivative of logΓ(.) and it is called ϕ(sk+1,sk ) e−∆E(skl,+n1)q(skl,n|skl,+n1) (15) l,n l,n ∝ q(sk+1sk ) l,n | l,n 4.4 Technical Details of the Algorithm where ∆E(sk+1) = E(sk+1) E(sk ) and E(sk ) = l,n l,n − l,n l,n We represent the proposed Adaptive Langevin Sampler W(sk ) + U(sk ). For any single pixel, U(s ) can be 1:L,n l,n l,n (ALS) algorithm in Appendix A. The symbol denotes derived from (3) and (6) as ←− analyticalupdate,thesymbol denotesupdatebyfind- 0 ←− D ing zero and the symbol denotes the update by ran- U(s )= 1+βl,d log 1+ φd(sl,n,αl,d) (16) dom sampling. The samplin∼g of the sources is done by the l,n 2 β δ Xd=1 (cid:20) l,d l,d (cid:21) Metropolis-Hastings scheme with Langevin proposal equa- tion. The random map produced by Langevin proposal is Theproposaldensityq(sk+1sk )isobtained,from(13), l,n | l,n appliedtoathresholdfunctiontokeeptheintensitiesofthe as maps in the physical margins. We have used the following τ sk+1sk + l,ng(sk ),τ2 (17) margins for CMB, synchrotron, dust and free-free, respec- N l,n | l,n 2 l,n l,n tively, [ 0.45,0.45], [0,0.5], [0,25] and [0,0.1]. They are in (cid:16) (cid:17) − The Metropolis-Hastings steps and the Langevin pro- antenna temperature (∆T)A in mK same as the maps. We posal equation are embedded into the main algorithm, as have determined these margins by using five patches from detailed in Appendix A. With this algorithm, we approach oursimulations. thesolution iteratively, avoidingtheinversion oftheconvo- lution matrix H and themixing matrix A. k 4.4.1 Initialization 4.3 Parameters of t-distribution By considering previous studies (Bonaldi et al. 2007; Ricciardi et al.2010),weassumethattheparametricmixing We find the mode estimates of the parameters of the t- matrixisknownwithanerror.Weinitializethemixingma- distribution using EM method. We can write the joint trixbyusingthespectralindicesobtainedin(Bonaldi et al. posterior of the parameters α , β and δ such that l,d l,d l,d 2007;Ricciardi et al.2010).Theparametricmodelisformed p(α ,β ,δ t ,Θ ) = p(t Θ)p(β ). In l,d l,d l,d| l,d −{αl,d,βl,d,δl,d} l,d| l,d sothatthecolumnsofsynchrotronanddustvaryaccording EM,ratherthanmaximizinglog p(t Θ)p(β ) ,wemax- l,d l,d to power laws only depending on one spectral index. Pre- { | } imize the following function iteratively vious experiments show that the error in spectral index of Θk+1=argmaxQ(Θ;Θk) (18) synchrotron changes from patch to patch and takes a max- Θ imum valueof about 1.72%. For spectral index of dust, the where superscript k representsthe iteration numberand maximum error is about 0.58%. We havefixedthe columns ofCMBandfree-freeastheyareknown.Weobtainrealistic Q(Θ;Θk)= log p(t Θ)p(β )) (19) h { l,d| l,d }iνl,d|tkl,d,Θk observations by mixing components with a mixing matrix which is formed by using the spectral indices 2.9 for syn- where p(ν tk ,Θk) is the posterior density of the hidden l,d| l,d chrotron and 1.8 for dust. In the reconstruction part, we variableν conditionedonparametersestimatedinthepre- l,d assumethatthespectralindicesareestimatedwithanerror vious step k and . represents the expectation hiνl,d|tkl,d,Θk of1.72% for synchrotronand0.58% for dust.So,weinitial- with respect toν tk ,Θk.Forsimplicity, hereafterweuse ize the mixing matrix values to maximum error case such l,d| l,d only the notation . torepresent this expectation. that the spectral indices are equal to 2.85 for synchrotron hi In the E (expectation) step of the EM algorithm, the and 1.7894 for dust. posterior expectation of ν is found as in Kayabol et al. To estimate the spectral indices, one can use the FD- l,d (2010) CCA (Fourier Domain Correlated Component Analysis) (Bedini & Salerno 2007) method, but it is not necessary to ν = N+βlk,d 1+ φd(skl,αkl,d) −1 (20) use this algorithm. Non-parametric mixing matrix estima- h l,di βk βk δk tionmethodscanbeused,suchasIndependentComponent l,d (cid:18) l,d l,d (cid:19) Joint Bayesian separation and restoration 5 Analysis(ICA)(Hyvarinen& Oja1997)orSpectralMatch- Gaussian shaped functions according to antenna apertures. ing ICA (SMICA) (Cardoso et al. 2002). Their standard deviations in pixels are given in Table 1. To initialize the component maps, we ignore the an- Weusethreedifferentperformancemeasuresdefinedin tennabeamsandapplytheinverseoftheinitialmixingma- thepixeldomaintoevaluatethesuccessoftheproposedal- trixtotherawobservationsdirectly.Ifwedenotetheinitial gorithm among the others. The Peak Signal-to-Interference mixing matrix A0, we initialize the maps with LS solution Ratio (PSIR ) in the pixeldomain is defined as pix as s0(n) = ((A0)TA0)−1(A0)Ty(n) where the vector y(n) √Nmax(s∗) containstheobservationintensitiesatnth pixel.LSsolution PSIR =20log l (25) pix RE is not a good solution since it does not take the noise and (cid:18) || || (cid:19) theresolutionsoftheobservationsintoconsideration,butit whereRE =s∗ ˆs istheReconstructionErrorbetweenthe provides a simple solution without any preprocessing inter- ground-truthsl∗−anldtheestimatedimageˆs .Weusethiser- l l vention. In this way, our algorithm starts with initial maps ror measure instead of the absolute difference between the which are some linear combination of theraw observations. ground-truth and the estimate, because we need a normal- Theinitialvaluesofαl,d canbecalculateddirectlyfromim- ized error to compare the error in large variations of the age differentials. Weinitialized the βl0,d =20 and found the intensitiesof thedifferentcomponents. Thelogarithm gives initial value of δl,d by equaling the expectation in (22) to a a good observation possibility for the values varying in a constant, in this study, we take it to be equal to 1.5. The large scale. Fig. 1, 2, 3, 4 and 5 show the estimated maps initial valueis found as δl0,d =1.5φd(s0l,α0l,d)/N. locatedatthecoordinatesof(0,0),(0,20),(0,40),(0,60)and (0,80). Wecompareourproposed methodALSwith theS+LS 4.4.2 Stopping Criterion andDB+LSsolutions.S+LSsolutionisobtainedbysmooth- We observe the normalized absolute difference ǫk = sk ing the observed data to the resolution of the lowest reso- sk−1/sk−1 between sequential values of s to dlecide| lth−e lution channel and then applying the inverse of the mixing l | | l | l matrix to the equal resolution maps. In DB+LS, we first convergence of the Markov Chain to an equilibrium. If the average ǫ¯k = 1 k ǫt 65e 2, we assume that thechain apply a de-blurring (DB) process to observation channels. l k t=1 l − For deblurring, we apply Wiener filter separately to each hasconvergedtotheequilibriumfors anddenotethispoint l T = k. Since wPe have L parallel chains for L sources, the channel using the known psf and the noise level. We find l the DB+LS solution by applying the inverse of the mixing endingpointoftheburn-inperiodofthewholeMonteCarlo matrix to the deblurred maps. chain isT =max T.Weignore thesamples before T .We s l l s For thepatch (0,0), we obtained a good reconstruction keep the iteration going until T that is the ending point e for CMB with theproposed method (Fig. 1). In themiddle of the post burn-in period simulation. In the experiments, of the map, the effect of the dust has been already seen. we have used 100 iterations after burn-in period, so T = e The effect of the dust also exists in the synchrotron map. T +100. s Thefree-freecomponentradiationmapinthepatches(0,0), At the end of the simulation the final estimates of the (0,20)and(0,40)cannotbeestimatedbyanymethod,since component maps are calculated as its intensity is very weak. For all the patches, the proposed 1 Te method reconstructs better the CMB and the foreground sˆl,n= T T skl,n (24) maps in the sense of PSIRpix. e s − kX=Ts Wealsoestimatetheerrorinthemaps.UsingMCsam- ples, we can find theuncertainty in the estimation. We call thisMCerrorandcalculateitsstandarddeviationforsingle pixelas 5 SIMULATION RESULTS 1 Inordertotestourideas,wehaveusedasetofrealisticsim- 1 Te 2 σ = (sˆ sk )2 (26) uoflamtiaopnssaonbdtatinoeodlsfdroemveltohpeedPlbaynctkheSPkylaMncokdWelo(rPkSinMg)G, arosuept MC Te−Ts kX=Ts l,n− l,n ! 2 (WG2) team as a fundamental part of the preparation where the numberT is the ending point of the simulation. e forthePlanckmission (Tauberet al.2010).Apartfrom the We use 100 iterations after convergence, so in this study CMB itself, the PSM contains state-of-the-art simulations T =T +100. e s of all the relevant Galactic and extragalactic astrophysical Weobtain anothererror measure byfittinga Gaussian components; for this work we use a simplified set of simu- to the posterior of the source image using the Laplace Ap- lations that contains CMB and Galactic (synchrotron,free- proximation (LA)method and calculating thestandard de- freeanddust)componentsonly,plusinstrumentalnoise.We viation, σ , of the approximated Gaussian. Table 2 lists LA have used simulations of the nine Planck frequencies. The the average standard deviations σ , σ and σ of the RE MC LA maincharacteristicsofthesimulationsarelistedinTable1. reconstruction, the Monte Carlo and the Laplace Approxi- We havetested ouralgorithm on fivedifferent 128x128 mation(LA)errors,respectively.Thereconstructionerroris patchesdistributedalongthecentralgalacticmeridian,and alwaysgreaterthantheestimatederrors.MCandLAerrors centered at galactic coordinates (00,00), (00,20), (00,40), arequiteclosetoeachother.TheminimumerrorsforCMB (00,60)and(00,80).Theactualsizeofthepatchesonthesky are found in thepatches (0,20) and (0,40). is 14.65◦ and thepixelsize is 6.87 arcmin. The mapsare in The plots in Fig. 6 compare the angular power spec- antennatemperature(∆T) inmK.Therelatednoiselevels trum, C defined as C = (ℓ+1)ℓC /2π, of the ground- A ℓ ℓ ℓ arepresentedinTable1.Wemodeltheblurringfunctionsas truth CMB maps and those obtained by S+LS, DB+LS 6 Kayabol et al. Groundtruth S+LS DB+LS ALS 0.2 00.2 0 0.2 CMB 0 −−00..42 −1 −00.2 −0.2 −0.6 −2 −0.8 −0.4 18.66 17.56 29.18 on 0.3 0.4 0.4 0.4 Synchrotr 00..12 000...123 −00.02.2 000...123 24.33 20.39 27.02 20 12 20 20 st 15 810 15 15 Du 10 6 10 10 5 4 5 5 2 35.77 102.4 111.3 x 10−4 4 0.4 0.1 −free 3 0.2 01.5 0.05 ee 2 0 Fr 1 0 −0.2 0 12.69 12.27 12.21 Figure 1. The estimated astrophysical maps at 100 GHz reference frequency from blurred and noisy observations with the S+LS, DB+LSandproposedALSmethods.Thelocationofthepatchis0◦ galacticlongitudeand0◦latitude.ThePSIRpixvaluesaredenoted undertheeachmap. Groundtruth S+LS DB+LS ALS 0.2 0.2 0 0.2 MB 0 0 −0.5 0 C −0.2 −1 −0.2 −0.2 24.94 25.49 41.2 n 0.04 0.08 0.04 o 0 otr 0.03 0.06 0.03 hr 0.04 −0.1 c Syn 0.02 0.02 −0.2 0.02 0 12.44 8.457 26.86 2.5 1.2 2.5 2.5 Dust 112.5 0001...468 112.5 112.5 0.5 0.2 0.5 0.5 41.37 86.31 93.5 x 10−4 0.05 0.6 0.025 e−free 24 0 00..24 00..00125 Fre −0.05 0 0.01 0 −0.1 0.005 15.84 13.06 18.18 Figure2.Theestimatedastrophysicalmapsat100GHzreferencefrequencyfromblurredandnoisyobservationswiththeS+LS,DB+LS andproposedALSmethods.Thelocationofthepatchis0◦ galacticlongitudeand20◦ latitude.ThePSIRpix valuesaredenotedunder theeachmap. Joint Bayesian separation and restoration 7 Groundtruth S+LS DB+LS ALS 0.2 0.2 0.2 0.2 B 0.1 CM 0 0 0 0 −0.2 −0.1 −0.2 −0.2 24.58 34.26 40.91 0.018 n 0.05 0.04 otro 0.016 0.04 0.02 0.02 Synchr 000...00011124 000...000123 −00.02 00..00115 4.858 4.381 18.11 0.03 0.02 0.025 0.03 0.02 ust 0.02 0.015 0.015 0.02 D 0.01 0.01 0.01 0.01 0.005 0.005 31.8 49.32 53.59 x 10−4 0.05 0.1 0.02 ee−free 510 0 00.05 0.015 Fr −0.05 −0.05 0.01 0 16.53 16.24 20.59 Figure3.Theestimatedastrophysicalmapsat100GHzreferencefrequencyfromblurredandnoisyobservationswiththeS+LS,DB+LS andproposedALSmethods.Thelocationofthepatchis0◦ galacticlongitudeand40◦ latitude.ThePSIRpix valuesaredenotedunder theeachmap. Groundtruth S+LS DB+LS ALS 0.2 0.2 0.2 0.2 0.1 MB 0 0 0 0 C −0.2 −0.1 −0.2 −0.2 23.95 33.09 40.61 0.02 0.03 n 0.06 0.06 o chrotr 0.015 0.04 00..0024 0.02 Syn 0.01 0.02 0 0.01 0 10.69 10.33 20.51 x 10−3 x 10−3 x 10−3 x 10−3 2 3 3 3 st 2.5 1.5 Du 2 2 1 2 1.5 1 0.5 1 26.76 27.16 35.4 0.08 0.1 0.08 ee−free 00..0046 000..015 0 00..0046 Fr 0.02 −0.05 0.02 −0.1 14.42 13.95 26.11 Figure4.Theestimatedastrophysicalmapsat100GHzreferencefrequencyfromblurredandnoisyobservationswiththeS+LS,DB+LS andproposedALSmethods.Thelocationofthepatchis0◦ galacticlongitudeand60◦ latitude.ThePSIRpix valuesaredenotedunder theeachmap. 8 Kayabol et al. Table 1.Channel frequencies, the standard deviations ofthe related Gaussian point spreadfunctions and thenoisestandarddeviations in(∆T)A [mK]. Channel frequencies[GHz] 30 44 70 100 143 217 353 545 857 psfstd .[pixels] 7.0069 4.8836 2.9726 2.1233 1.5075 1.0617 1.0617 1.0617 1.0617 Noisestd (∆T)A [mK] 0.0259 0.0248 0.0233 0.0074 0.0038 0.0032 0.0023 0.0019 0.0009 Groundtruth S+LS DB+LS ALS 0.2 0.2 0.2 0.2 B M 0 0 0 0 C −0.2 −0.2 −0.2 −0.2 24.79 34.5 41.22 0.06 0.06 hrotron 000...000111246 0.04 00..0024 00..00125 nc 0.01 0.02 0 0.01 Sy 0.008 −0.02 0.006 0 0.005 10.64 10.09 20.85 x 10−3 x 10−3 x 10−3 x 10−3 6 4 6 3 3 Dust 4 2 2 4 2 1 2 1 32.16 33.64 41.69 0.08 0.08 0.1 0.1 ee−free 00..0046 00.05 00.05 00..0046 Fr 0.02 −0.05 −0.05 0.02 15.89 15.56 26.13 Figure5.Theestimatedastrophysicalmapsat100GHzreferencefrequencyfromblurredandnoisyobservationswiththeS+LS,DB+LS andproposedALSmethods.Thelocationofthepatchis0◦ galacticlongitudeand80◦ latitude.ThePSIRpix valuesaredenotedunder theeachmap. and the proposed ALS methods. To plot C , we sample it petitor methods. The algorithm works quite well at high ℓ by taking √N/2+1 samples in the interval [0,ℓ ] where latitudes. If we approach the galactic plane, the estimation max ℓ = 180/14.65√N. All the patches, the power spectra results get worse. Especially at the galactic plane, we have max found by the proposed ALS method fit the ground-truth obtainedtheworstresults,althoughwehaveusedadifferent spectra more tightly than the others. The S+LS method initialization strategy. gives bad results in the high frequency regions of the spec- Our new goal is the application of the proposed algo- trums because of smoothing. The DB+LS method causes rithmtowhole-skymaps.Toavoidthedifficultiesinherentin anattenuationinthelowfrequencyregions.Therootmean thisproblem,weplan tousethe”nestednumbering”struc- square error between the groundtruth angular power spec- tureprovidedbytheHEALPix(Gorski et al.2005)package. trum of CMB and those obtained by S+LS, DB+LS and In this format, we can reach theindexes of theeight neigh- proposed ALS methods are presented in Table 3. The pro- borsofeachpixelonthesphere.Tocalculatethepixeldiffer- posedmethodprovideoneorderofthemagnitudebetterfit ences, we will implement a gradient calculation method on than theothers. thespherebytakingthenon-homogeneousspatialdistances between the pixels on thesphere into consideration. 6 CONCLUSIONS ACKNOWLEDGMENTS WehaveintroducedaBayesianjointseparationandestima- tion method for astrophysical images. The method is based The authors would like to thank Anna Bonaldi,(INAF, onaMonteCarlotechniqueandgivesbetterreconstruction Padova,Italy),BulentSankur,(BogaziciUniversity,Turkey) in the pixel domain and frequency domain than two com- and Luigi Bedini (ISTI, CNR, Italy) for valuable discus- Joint Bayesian separation and restoration 9 x 104 CMB CMB Groundtruth 12000 Groundtruth 2 S+LS S+LS DB+LS DB+LS ALS 10000 ALS 1.5 8000 π π 2 2 C/ C/ +1)l 1 +1)l 6000 (l (l 4000 0.5 2000 0 0 0 500 1000 1500 0 500 1000 1500 l l (a) Patch(0◦,0◦) (b) Patch(0◦,20◦) CMB CMB 10000 Groundtruth Groundtruth S+LS 10000 S+LS 9000 DB+LS DB+LS ALS 9000 ALS 8000 8000 7000 7000 πC/2 6000 πC/2 6000 +1)l 5000 +1)l 5000 (l (l 4000 4000 3000 3000 2000 2000 1000 1000 0 0 0 500 1000 1500 0 500 1000 1500 l l (c) Patch(0◦,40◦) (d) Patch(0◦,60◦) CMB 10000 Groundtruth 9000 S+LS DB+LS 8000 ALS 7000 6000 π 2 1)lC/ 5000 + (l 4000 3000 2000 1000 0 0 500 1000 1500 l (e) Patch(0◦,80◦) Figure 6. Comparison of the standard power spectrum of the ground-truth maps located at a) (0◦,0◦), b) (0◦,20◦), c) (0◦,40◦), d) (0◦,60◦)ande)(0◦,80◦)withthoseobtainedbyS+LS,DB+LSandproposedALSmethods. 10 Kayabol et al. is partially supported by CNR-CSIC bilateral project no: Table 2. Average standard deviations of Reconstruction Error (RE),MonteCarlo(MC)uncertaintyandLaplaceApproximation 2008IT0059. (LA)uncertainty. (0◦,0◦) CMB Synchrotron Dust Free-Free REFERENCES σRE 30.49e-3 13.82e-3 15.04e-3 16.32e-3 σMC 2.38e-3 0.42e-3 0.08e-3 1.50e-3 AnthoineS., 2005, PhD Thesis, Princeton University σLA 3.38e-3 0.14e-3 0.06e-3 1.28e-3 BeckerS.,LeCunY.,1989,Proc.ofthe1988Connectionist (0◦,20◦) Models Summer School, 29 Bedini L., Herranz D., Salerno E., Baccigalupi C., Ku- CMB Synchrotron Dust Free-Free ruog˘lu E.E., Tonazzini A., 2005, EURASIP Journal on σRE 14.83e-3 2.26e-3 1.03e-3 12.80e-3 Applied Signal Processing, 15, 2400 σMC 2.81e-3 1.04e-3 0.03e-3 2.10e-3 Bedini L., Salerno E., 2007, Lecture Notes in Artificial In- σLA 3.37e-3 0.92e-3 0.06e-3 1.46e-3 telligence, 4694, 9 Belouchrani, A.,A.-Meraim, K.,Cardoso, J.-F., Moulines, (0◦,40◦) E., 1997, IEEE Trans. Signal Process., 45, 434 CMB Synchrotron Dust Free-Free Bonaldi A., Ricciardi S., Leach S., Stivoli F., Baccigalupi C., De Zotti G., 2007, MNRAS,382, 1791 σRE 14.22e-3 2.47e-3 0.17e-3 12.25e-3 Cardoso J.-F., Snoussi H., Delabrouille J., Patanchon σMC 2.88e-3 1.28e-3 0.03e-3 1.96e-3 G., 2002, European. Conf. on Signal Processing, EU- σLA 3.37e-3 1.12e-3 0.02e-3 1.48e-3 SIPCO’02, 561 (0◦,60◦) Castella M., Pesquet J.-C., 2004, LNCS,ICA, 3195, 922 Chantas G., Galatsanos N., Likas A., Saunders M., 2008,, CMB Synchrotron Dust Free-Free IEEE Trans. Image Process., 17, 1795 σRE 10.72e-3 2.85e-3 0.06e-3 8.27e-3 Eriksen H.K., Jewel J.B., Dickinson C., Banday A.J., σMC 2.72e-3 1.34e-3 0.02e-3 2.05e-3 Gorski K.M., Lawrence C.R., 2008, ApJ, 676, 10 σLA 3.37e-3 1.05e-3 0.02e-3 1.53e-3 Gorski, K.M., Hivon, E., Banday, A.J., Wandelt, B.D., (0◦,80◦) Hansen,K.M., Reinecke, F., Bartelmann, M., 2005, ApJ, 622, 759 CMB Synchrotron Dust Free-Free Hastings W.K., 1970, Biometrika, 57, 97 Higdon D., 1994, PhD Thesis, University of Washingthon σRE 10.29e-3 3.21e-3 0.06e-3 7.40e-3 σMC 2.83e-3 1.43e-3 0.02e-3 2.11e-3 Hyvarinen,A.,Oja,E.,1997,NeuralComputation,9,1483 σLA 3.37e-3 1.11e-3 0.02e-3 1.54e-3 KayabolK.,KuruogluE.E.,SankurB.,2009,IEEETrans. Image Process., 18, 982 KayabolK., KuruogluE.E., SanzJ.L., SankurB., Salerno Table3.Rootmeansquareerrorbetweenthegroundtruthstan- E., Herranz D., 2010, IEEE Trans. Image Process. dard power spectrum of CMB and those obtained by S+LS, Liu C., Rubin D.B., 1995, Statistica Sinica, 5, 19 DB+LSandproposedALSmethodsatpatches(0◦,0◦),(0◦,20◦), Maino D., Farusi A., Baccigalupi C., Perrotta F., Banday (0◦,40◦),(0◦,60◦)and(0◦,80◦). A.J., Bedini L., Burigana C., De Zotti G., G´orski K.M., Salerno E., 2002, MNRAS,334, 53 S+LS DB+LS ALS Prudyus I., Voloshynovskiy S., Synyavskyy A., 2001, Int. Conf.onTelecomm. inModernSatell., CableandBroad- (0◦,0◦) 2.1798e+3 7.4638e+3 0.1894e+3 cas. TELSIKS’01, 583 (0◦,20◦) 2.1867e+3 0.5742e+3 0.0605e+3 (0◦,40◦) 2.2902e+3 1.5089e+3 0.0436e+3 RicciardiS.,BonaldiA.,NatoliP.,PolentaG.,Baccigalupi (0◦,60◦) 2.2261e+3 1.5380e+3 0.0462e+3 C., Salerno E., KayabolK.,Bedini L., DeZottiG., 2010, (0◦,80◦) 2.1598e+3 1.4089e+3 0.0513e+3 MNRAS, Shwartz S., Schechner Y.Y., Zibulevsky M., 2008, Neuro- computing, 71, 2164 sions.ThesimulatedsourcemapsaretakenfromthePlanck Tauber, J. A., Mandolesi, N., Puget, J.-L., Banos, T., Sky Model, a set of maps and tools for generating realistic Bersanelli, M., Bouchet, F. R., Butler, R. C., Charra, J., Planck simulations made available thanks to the efforts of Crone, G., Dodsworth, J., and et al., 2010, A&A,520, 1 the Planck Working Group 2 (WG2) team. Some of the re- Tonazzini A.,Gerace I., 2005, European Conf. Signal Pro- sults in this paper have been derived using the HEALPix cess., EUSIPCO Gorski et al. (2005) package. Tonazzini A., Gerace I., Martinelli F., 2010, IEEE Trans. KorayKayabolundertookthisworkwiththesupportof Image Process., 19, 912 the”ICTPProgrammeforTrainingandResearchinItalian TzikasD.,LikasA.,GalatsanosN.,2009,IEEETrans.Im- Laboratories, Trieste, Italy, through a specific operational age Process., 18, 753 agreement with CNR-ISTI, Italy. Partial support has also Wilson S.,KuruogluE.E., SalernoE., 2008, IEEEJournal beengivenbytheItalianSpaceAgency(ASI),underproject on Selected Areas in Signal Processing, 2, 685 COFIS(CosmologyandFundamentalPhysics).Theproject