“Internal Linear Combination” method for the separation of CMB from Galactic foregrounds in the harmonic domain
Key Words.:
Methods: data analysis – Methods: statistical – Cosmology: cosmic microwave backgroundForeground contamination is the fundamental hindrance to the cosmic microwave background (CMB) signals and its separation from it represents a fundamental question in Cosmology. One of the most popular algorithm used to disentangle foregrounds from the CMB signals is the “internal linear combination” method (ILC). In its original version, this technique is applied directly to the observed maps. In recent literature, however, it is suggested that in the harmonic (Fourier) domain it is possible to obtain better results since a separation can be attempted where the various Fourier frequencies are given different weights. This is seen as a useful characteristic in the case of noisy data. Here, we argue that the benefits of using such an approach are overestimated. Better results can be obtained if a classic procedure is adopted where data are filtered before the separation is carried out.
1 Introduction
The experimental progresses in the detection of cosmological emissions require a parallel development of data analysis techniques in order to extract the maximum physical information from data. In particular, different emission mechanisms are characterized by markedly distinct underlying physical processes. Data analysis often requires the component separation in order to study the individual characteristics. To achieve such a goal, a link between the branch of signal processing science which characterizes and separates different signals and astrophysics is yet well established, and in many cases, modern signal processing techniques have been imported and applied in an astrophysical context.
In this work, we consider one of the most widely used approaches for the separation of different emissions, say the internal linear combination method (ILC), which was adopted, for instance, in the reduction of the data from the Wilkinson Microwave Anisotropy Probe (WMAP) satellite for CMB observations (Bennett et al. 2003). Among the separation techniques, ILC requires the smallest number of a priori assumptions. Here, the available data are assumed in the form of maps, taken at different frequencies, containing pixels each. More precisely, if provides the value of the th pixel for a map obtained at channel “” ^{1}^{1}1In the present work, indexes pixels in the classic spatial domain. However, the same formalism applies if other domains are considered, for example, the Fourier one., our starting model is:
(1) 
where , and are the contributions due to the CMB, the diffuse Galactic foreground and the experimental noise, respectively. Although not necessary for later arguments, it is assumed that all of these contributions are representable by means of stationary random fields. Moreover, without loss of generality, for ease of notation the random fields are supposed as the realization of zeromean spatial processes. In the present work the contribution of nondiffuse components (e.g., due to SZ cluster, pointsources, …) are not considered and they are assumed to have been removed through other methodologies.
The idea behind ILC is simple. The main assumption is that model (1) can be written as
(2) 
i.e. the template of the CMB component is independent of the observing channel. A way to exploit this fact is to average images giving a specific weight to each of them so as to minimize the impact of the foreground and noise (Bennett et al. 2003). This means to look for a solution of type
(3) 
If the constraint is imposed, Eq. (3) becomes
(4) 
Now, from this equation it is clear that, for a given pixel “”, the only variable terms are in the sum. Hence, under the assumption of independence of from and , the weights have to minimize the variance of , i.e.
(5) 
where is the expected variance of . If denotes a row vector such as and the matrix is defined as
(6) 
then Eq. (2) becomes
(7) 
In this case, the weights are given by (Eriksen et al. 2004)
(8) 
where is the crosscovariance matrix of the random processes that generate , i.e.
(9) 
and is a column vector of all ones. Here, denotes the expectation operator. Hence, the ILC estimator takes the form
(10)  
(11) 
with and the scalar quantity given by
(12) 
In practical applications, matrix is unknown and has to be estimated from the data. Typically, this is done by means of the estimator
(13) 
A common assumption in CMB observations is that is given by the linear mixture of the contribution of physical processes (e.g., freefree, dust reradiation, …)
(14) 
with constant coefficients. In practice, it is assumed that for the th physical process a template exists independent of the specific channel “”. Although rather strong, it is not unrealistic to assume that this condition is satisfied when small enough patches of the sky are considered. Inserting Eq. (14) into Eq. (7) one obtains
(15) 
with
(16) 
and
(17) 
Here, matrix is assumed to be of full rank.
2 ILC in the spatial domain
In a recent paper, Vio & Andreani (2008) have stressed various problems concerning ILC that in literature have been underevaluated if not undetected. For example, in Eq. (7) the term is often considered as a single noise component (e.g., see Eriksen et al. 2004; Hinshaw et al. 2007; Kim et al. 2007, 2008). In this way the problem is apparently simplified since it is reduced to the separation of two components only. No a priori information on this “global” noise is required. However, this approach can lead to wrong conclusions. For example, since all the components in the mixtures are assumed to be zeromean, from Eq. (4) one could conclude that
(18) 
i.e. the ILC estimator is unbiased ^{2}^{2}2The expression indicates conditional expectation of given . . This is not correct: the claim that is unbiased requires one to prove that
(19) 
The reason is that is a fixed realization of a random process. There is no way to obtain another one. Even if observed many times (under the same experimental conditions) the foreground components (for instance the Galaxy) will always appear the same. Only the noise component will change. This has important consequences. In fact, in Vio & Andreani (2008) it is proved that, even under model (15), ILC can provide satisfactory results only under rather restrictive conditions:

The number of observing frequencies has to be larger than the number of components ;

has to be uncorrelated with ;

Data have to be noisefree, i.e. .
The violation of any of these points has as consequence the introduction of a bias in the ILC solution that can be severe. In literature it appears that only the second point has been well realized. (e.g. see Delabrouille & Cardoso 2007). The explanation of the bias when is technical and we remand to Vio & Andreani (2008). When noise is present, under the condition of uncorrelated with , from model (15) it is
(20) 
with
(21) 
Hence, a bias derives from the fact that is a matrix with strictly positive entries and then .
In absence of “a priori information”, the problems connected to the first two points above have no solution. The only possibility is to plan the experiments in such a way to avoid them (e.g. observations at high Galactic latitudes, a number of observing frequencies sufficiently large …). On the other hand, noise is an unavoidable question that, however, can be handled with hope of satisfactory results. In this respect, some authors (e.g. see Tegmark et al. 2003; Kim et al. 2007, 2008) suggest that an effective way is to use ILC in the harmonic (Fourier) domain. In practice, this means to apply ILC to that is the matrix whose th row contains the twodimensional Fourier transform of . Following Kim et al. (2007), we will indicate this version of ILC as “harmonic internal linear combination” (HILC).
3 ILC in the Fourier domain (HILC)
The starting consideration is that all arguments presented in the previous sections hold independently of the fact that one is working in the ordinary spatial domain or in the Fourier domain. Indeed, it is not difficult to see that the weights as given by Eq. (8) can be equivalently computed by means of
(22) 
In this respect, it is sufficient to take into account that
(23) 
with , “” the Kronecker product, the the Fourier matrix that is a complex, symmetric and unitary matrix whose elements are given by
(24) 
. Hence
(25) 
since where is the complex conjugate transpose of .
It is well known that measurement noise , typically the realization of wideband stochastic processes, tends to uniformly spread in the Fourier domain. On the other hand, the contribution of signals is typically concentrated in correspondence to the lowest Fourier frequencies. For example, this is visible in Fig. 1 where the power spectra of three simulated mixtures are shown (see Fig. 3). Here, nonastronomical signals have been used, but the same holds also for astronomical ones. This implies that, for the indices “” corresponding the lowest frequencies, the contribution of noise to is quite small. Hence, the basic idea here is to partition the frequency domain in a number of subsets and to apply separately the ILC filter to each of them. At the end of this procedure an estimate is obtained where . Finally, can be recovered by the Fourier inversion of .
The effectiveness of such an approach seems supported by the simple example presented in Figs. 24. In particular, Fig. 2 shows the original images , , as well the mixing matrix that, through model (15), are used to create the observed images shown in Fig. 3. As reference for later results, the bottomright panel in the same figure shows the estimate obtained when ILC is applied to these images. For this case, the resulting root mean square (rms) of the residuals is about . The top panels in Fig. 4 show what happens when the images are added a whitenoise with a signaltonoise (SNR) ^{3}^{3}3Here, the quantity SNR is defined as the ratio between the standard deviation of the signal with the standard deviation of the noise. set to . Here, both ILC and HILC are used. In the case of HILC, the frequency domain is partitioned in two regions, as shown in the bottomleft panel of Fig. 1, one containing the low frequencies and the other containing the high frequencies. The weights are calculated independently for each of them. It is evident that HILC outperforms ILC. This conclusion is confirmed by the rms of the its residuals that is and , respectively. The same indication is provided by the bottom panels that show the result obtained when the weights computed for the noisy images are applied to the noisefree ones. This operation gives an idea of the bias introduced in the solution by the noise. In this case, the rms of the residuals become and , respectively.
4 Discussion and conclusions
Reader should not be confused by the fact that at first sight HILC appears an approach potentially superior to ILC. It is indeed true that HILC performs a more effective separation of the signal from the noise. However, if the noise affects most the high frequency part of the signal, why then do not simply filter out this one? Indeed, the left panel of Fig. 5 shows the estimate obtained when, before the application of ILC, the images are filtered with the ideal Fourier circular lowpass filter ^{4}^{4}4In the Fourier domain, an ideal filter has the form if is a frequency of interest, otherwise. that has the structure identical to that shown in Fig. 1. As before, the right panel of this figure shows the solution obtained when the resulting weights are applied to the noisefree version of the data. The rms of the residuals are and , respectively. The comparison with the values of and , previously determined for HILC, indicates an improvement of the quality of the result. The indication that comes out from this simple experiment is that filtering noise is a more effective operation than its separation from the signal of interest. It is true that the version of HILC here used is rather rough and more sophisticated ones are available (Kim et al. 2007, 2008). However, the same holds also for the filtering operation that has been coupled with ILC. At this point a question arises: which is the benefit in using a nonstandard approach as HILC with respect to a classic approach where filtering is coupled with ILC? This in not an academic question. Indeed, the use of nonstandard tools implies the renunciation of a huge body of experience gained in years of application in practical problems also in disciplines different from Astronomy. In other words, the choice of nonstandard tools is indicated only in situations of real and sensible improvements of the results. New techniques that do not fulfill this requirement should be introduced with care: they prevent the comparison with the results obtained in other works and may lead people to use not well tested methodologies ending up in not reliable results.
References
 Bennett et al. (2003) Bennett, C.L., et al. 2003, ApJS, 148, 97
 Delabrouille & Cardoso (2007) Delabrouille, J., & Cardoso, J.F. 2007, arXiv:astroph/0702198v2
 Eriksen et al. (2004) Eriksen, H.K., Banday, A.J., Górski, K.M., & Lilje, P.B. 2004, ApJ, 612, 633
 Hinshaw et al. (2007) Hinshaw, G., et al. 2007, ApJS, 170, 288
 Kim et al. (2007) Kim, j., Naselsky, P., & Christensen, P.R. 2007, Phys. Rev. D, 375, 625
 Kim et al. (2008) Kim, j., Naselsky, P., & Christensen, P.R. 2008, arXiv:0810.4008
 Tegmark et al. (2003) Tegmark, M., de OliveiraCosta, A., & Hamilton, A.J. 2003, Phys. Rev. D., 68, 123523
 Vio & Andreani (2008) Vio, R., & Andreani, P. 2008, A&A 487, 775