Table Of ContentMon.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
