Astronomy & Astrophysics manuscript no. (DOI: will be inserted by hand later)
17th January 2014
Weak Lensing Mass Reconstruction using Wavelets Jean-Luc Starck, Sandrine Pires, and Alexandre R´efr´egier DAPNIA/SEDI-SAP, Service d’Astrophysique, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France 17th January 2014 Abstract This paper presents a new method for the reconstruction of weak lensing mass maps. It uses the multiscale entropy concept, which is based on wavelets, and the False Discovery Rate (FDR) which allows us to derive robust detection levels in wavelet space. We show that this new restoration approach outperforms several standard techniques currently used for weak shear mass reconstruction. This method can also be used to separate E and B modes in the shear field, and thus test for the presence of residual systematic effects. We concentrate on large blind cosmic shear surveys, and illustrate our results using simulated shear maps derived from N-Body ΛCDM simulations (Vale and White, 2003) with added noise corresponding to both ground-based and space-based observations.
Key words. Cosmology : Weak Lensing, Methods : Data Analysis
1. Introduction Weak Gravitational Lensing provides a unique method to map directly the distribution of dark matter in the universe (Bartelmann and Schneider, 1999; Mellier, 1999; Van Waerbeke et al., 2001; Mellier, 2002; Refregier, 2003). This method is based on the weak distortions that lensing induces in the images of background galaxies as light travels through intervening structures. This method is now widely used to map the mass of clusters and superclusters of galaxies and to measure the statistics of the cosmic shear field on large scales. Ongoing efforts are made to improve of cosmic shear on existing telescopes and ments dedicated to survey cosmic shear Several methods are used to derive the Send offprint requests to: [email protected]
the detection future instruare planned. lensing shear
from the shapes of background galaxies. But the shear measurements obtained are always noisy, and when it is converted into a map of the projected mass κ, the result is dominated by the noise. Several methods have been devised to reconstruct the projected mass distribution from the observed shear field. The first non-parametric mass reconstruction was proposed by Kaiser and Squires (1993) and further improved by Bartelmann (1995); Kaiser (1995); Schneider and Seitz (1995); Squires and Kaiser (1996). These methods are based on linear inversion methods based on smoothing with a fixed kernel. Non-linear reconstruction methods were proposed using a maximum likelihood approach (Squires and Kaiser, 1996; Bartelmann et al., 1996; Seitz et al., 1998) or using the maximum entropy method (Bridle et al., 1998; Marshall et al., 2002). In this paper, we describe a method for weak lensing mass reconstruction based on a wavelet decomposition. We use an iterative filtering method with a multiscale entropy regularisation to filter the noise. We discuss how this decomposition and regularisation functional
Weak Lensing Mass Reconstruction
is particularly well adapted to this problem. In the process, we identify significant wavelet coefficients using the False Discovery Rate method (Miller et al., 2001; Hopkins et al., 2002) and show how this is superior to the standard nσ thresholding. The FDR method adapts its threshold to the features of the data. We concentrate on large blind cosmic shear surveys and use the raytracing simulations of Vale and White (2003) to test our results. We compare the performance of our method to Gaussian and Wiener filtering for the reconstruction of the mass field in these simulations. We consider conditions similar to both ground-based and space-based cosmic shear surveys. We also discuss how our method differs from other methods based on the maximum entropy prior. In section 2, we present the weak shear mass reconstruction problem. The earlier methods which have been proposed to reconstruct the mass map are described in section 3. Section 4 presents the Multiscale Entropy method and explains why it is a good alternative to standard methods. We also propose a modification of the Multiscale Entropy, we use the False Discovery Rate (FDR) method for detecting the significant wavelet coefficients. A set of experiments designed to test our method are described section 5. Our conclusions are summarized in section 6.
2. Weak lensing mass reconstruction 2.1. Weak lensing In weak lensing surveys, the shear γi (θ) with i = 1, 2 is derived from the shapes of galaxies at positions θ in the image. The shear field γi (θ) can be written in terms of the lensing potential ψ(θ) as (see eg. Bartelmann and Schneider (1999)) 1 2 γ1 = ∂1 − ∂22 ψ 2 γ2 = ∂1 ∂2 ψ, (1) where the partial derivatives ∂i are with respect to θi . The convergence κ(θ) can also be expressed in terms of the lensing potential as 1 2 ∂ + ∂22 ψ (2) κ= 2 1 and is related to the surface density Σ(θ) projected along the line of sight by κ(θ) =
where the critical surface density is given by Σcrit =
c2 D s 4πG Dl Dls
and G is Newton’s constant, c is the speed of light and Ds , Dl and Dls are the angular-diameter distances between the observer and the galaxies, the observer and the lens, and the lens and the galaxies. In practice, the galaxies are not at a fixed redshift, and the expression for κ is an average of the redshift of the galaxies (see eg. Bartelmann (1995)). The lensing effect is said to be weak or strong if κ ≪ 1 or κ & 1, respectively. The left panel of Fig. 1 shows a simulated convergence map derived from ray-tracing through N-body cosmological simulations performed by Vale and White (2003). The cosmological model is taken to be a concordance ΛCDM model with parameters ΩM = 0.3, ΩΛ = 0.7, h = 0.7 and σ8 = 0.8. The simulation contains 5123 particles with a box size of 300h−1 Mpc. The resulting convergence map covers 2 × 2 degrees with 1024 × 1024 pixels and a assume a galaxy redshift of 1. The overdensities correspond to the haloes of groups and clusters of galaxies. The rms value of κ binned in 0.12 arcmin pixels is σκ = 0.023. The typical values of κ are thus of the order of a few percent, apart from the core of massive halos (see figure 1). The weak lensing condition therefore holds in most regions of the sky and will be assumed throughout this paper.
2.2. Mass inversion The weak lensing mass inversion problem consists of reconstructing the projected (normalized) mass distribution κ(θ) from the measured shear field γi (θ) by inverting equations (1) and (2). (Magnification information can also be used to improve the reconstruction [see Bridle et al. (1998)], but is typically more noisy than the shear measurements and has not been considered in this paper). For this purpose, we take the Fourier transform of these equations and obtain γˆi = Pˆi κ ˆ,
i = 1, 2
where the hat symbol denotes Fourier transforms and we have defined k 2 ≡ k12 + k22 and k2 − k2 Pˆ1 (k) = 1 2 2 k 2k k2 1 Pˆ2 (k) = , k2
Figure 1. Left: simulated convergence map from (Vale and White, 2003) for a ΛCDM model. The region shown is 2 × 2 square degree. Right: Shear map superimposed on the convergence map , and right shear map. The size and direction of each line gives the amplitude and position angle of the shear at this location on the sky. with Pˆ1 (k1 , k2 ) ≡ 0 when k12 = k22 , and Pˆ2 (k1 , k2 ) ≡ 0 when k1 = 0 or k2 = 0. The shear map γi can be calculated from the convergence map κ using these expressions. The right panel of Fig. 1, shows the shear field associated with the simulated convergence field. As is customary, the direction and size of the line segment represent the orientation and amplitude of the shear. The rms shear in the 0.12 amin pixels 1 of the resulting map is σγ = hγ12 + γ22 i 2 ≃ 0.023. Note that to recover κ from γ1 (resp. γ2 ), there is a degeneracy when k12 = k22 (resp. when k1 = 0 or k2 = 0). To recover κ from both γ1 and γ2 , there is a degeneracy only when k1 = k2 = 0. Therefore, the mean value of κ cannot be recovered from the shear maps. This is a special instance of the well known mass-sheet degeneracy in the weak lensing reconstruction if only shear information is available (see eg.Bartelmann (1995) for a discussion). In practice, the observed shear γi is obtained by averaging over a finite number of galaxies and is therefore noisy. The relations between the observed data γ1b , γ2b binned in pixels of area A and the true mass map κ are given by: γib = Pi ∗ κ + Ni
where N1 and N2 are noise contributions with zero mean p and standard deviation σn ≃ σǫ / Ng , where Ng = ng A is the average number of galaxies in a pixel and ng is the average number of galaxies per arcmin2 . The rms shear dispersion per galaxy σǫ arises both from measurement errors and the intrinsic shape dispersion of galaxies. In this analysis, we will assume σǫ ≃ 0.3 as is approximately found for ground-based and space-based weak lensing surveys. Typical values for the surface density of usable galaxies for weak lensing are – ng = 20 gal/arcmin2 for ground-based surveys. – ng = 100 gal/arcmin2 for space-based surveys. From the central limit theorem, this means that for pixels with A & 1 amin2 , the noise Ni is, to a good approximation, Gaussian in both cases and is uncorrelated (see Marshall et al. (2002) for a direct treatment of individual galaxy shears using the MEM method).
2.3. The Inverse Filter: E and B mode We can easily derive an estimation of the mass map by inverse filtering by noticing that 2 2 Pˆ1 + Pˆ2 = 1.
ˆ˜ l(E) of the convergence κ The least square estimator κ ˆ in the Fourier domain is: (E) κ ˜ˆ l = Pˆ1 γˆ1b + Pˆ2 γˆ2b
The relation between this estimator and the true mass ˆ , where N ˆ = Pˆ1 Nˆ1 + Pˆ2 Nˆ2 . ˆ˜ l(E) = κ map is κ ˆ+N Just as any vector field, the shear field γi (θ) can be decomposed into a gradient, or electric (E), component, and a curl, or magnetic (B), component. Because the weak lensing arises from a scalar potential (the Newtonian potential), it can be shown that weak lensing only produces E-modes. On the other hand, residual systematics arising from imperfect correction of the instrumental PSF or telescope aberrations, generally generates both E and B modes. The presence of B-modes is thus used to test for the presence of residual systematic effects in current weak lensing surveys. The decomposition of the shear field into each of these components can be easily performed by noticing that a pure E-mode can be transformed into a pure B mode by a rotation of the shear by 45◦ : γ1 → −γ2 , γ2 → γ1 . As a result, we can form the following estimator for the B-mode “convergence” field ˆ˜ (B) κ = Pˆ2 γˆ1b − Pˆ1 ∗ γˆ2b , l
and check that it is consistent with zero in the absence of systematics.
Figure 2. In the previous region of 2 × 2 square de(E) grees, noisy mass map κl for the same simulation with ng = 100 gal/arcmin2, corresponding to space-based observations. Even in this case, the unfiltered mass map is dominated by noise.
Weak Lensing Mass Reconstruction
As follows from equation (8), the noise N (E) and N (B) (E) (B) in κ ˜l and κ ˜l is still Gaussian and uncorrelated. The in(E) verse filtering does not amplify the noise, but κ ˜ l and (B) κ ˜ l may be dominated by the noise if N (E) and N (B) are large, which is the case in practice. Fig. 2 shows the reconstructed mass map using equation 9 when a realistic Gaussian noise has been added to the shear maps plotted in Fig. 1 right. As expected, it is dominated by noise. This has motivated the development of different methods in the past which we describe below.
3.2. Maximum Entropy Method The Maximum Entropy Method (MEM) is well-known and widely used in image analysis in astronomy (see Bridle et al. (1998); Starck et al. (2001); Marshall et al. (2002); Starck and Murtagh (2002) for a full description). It considers both the data and the solution as probability density functions and find the solution using a Bayesian approach and adding a prior (the entropy) on the solution. Several definitions of entropy exists. The most common is the definition proposed in Gull and Skilling (1991): Hg (κ) =
3. Earlier Mass Inversion Methods
κ(x, y) − m(x, y) − κ(x, y) ln
κ(x, y) m(x, y)
3.1. Linear Filtering The standard method (Kaiser and Squires, 1993) con(E) sists in convolving the noisy mass map κ ˜ l with a Gaussian window G with standard deviation σG : (E)
κ ˜G = G ∗ κ ˜l
= G ∗ P1 ∗ γ1b + G ∗ P2 ∗ γ2b
The quality of the resulting estimation depends strongly on the value of σG . Fig. 3 shows the variation of the error between the original mass map κ shown in Fig.1 and the (E) filtered mass map κ ˜ G . For this simulation, the optimal value of σG lies between 5 and 10 pixels (1 pixel = 0.12 arcmin) for space observations (i.e. ng = 100 gal/amin2) and lies between 20 and 25 pixels for ground observations ( i.e. ng = 20 gal/amin2). An alternative to Gaussian filtering is the Wiener filtering, obtained by assigning the following weight to each k-mode w(k) =
|ˆ κ(k)|2 ˆ (k)|2 |ˆ κ(k)|2 + |N
Where |ˆ κ(k)|2 is a model of the true convergence power spectrum and is in practice derived from the data. Wiener filtering is known to be optimal when both the signal and the noise are a realization of a Gaussian Random Field. As can be seen from Fig.1, this assumption is not valid for weak lensing mass maps which display non-Gaussian features such as galaxy clusters, groups and filaments. Even in this case, Wiener Filtering nevertheless leads to reasonable results, generally better than the simple Gaussian filtering.
where m is a model, chosen typically to be a sky background. Hg has a global maximum at κ = m. MEM does not allow negative values in the solution, which is unnatural for wide field weak shear data or the CMB data, where we measure fluctuations around zero. To overcome this, it has been proposed to replace Hg (Maisinger et al., 2004) by: H+/− (κ) =
ψ(x, y) − 2m −
−κ(x, y) ln
ψ(x, y) + κ(x, y) 2m
p where ψ(x, y) = κ2 (x, y) + 4m2 . Here m does not play the same role. It is a constant fixed to the expected signal rms. More generally MEM method presents many drawbacks (Narayan and Nityananda, 1986; Starck et al., 2001) and various refinements of MEM have been proposed over the years (Weir, 1992; Bontekoe et al., 1994; Pantin and Starck, 1996; Starck et al., 2001). The last developments have lead to the so called Multiscale Entropy (Pantin and Starck, 1996; Starck et al., 2001; Maisinger et al., 2004) which is based on an undecimated isotropic wavelet transform (` a trous algorithm) (Starck et al., 1998). It has been shown that the main MEM drawbacks (model dependent solution, oversmoothing of compact objects, . . . ) disappear in the wavelet framework. A full discussion and comparison between different restoration methods can be found in Starck et al. (2002).
Figure 3. Reconstruction error as a function of the kernel size σG (in 0.12 amin pixels) for the Gaussian smoothing method, with ng = 20 gal/amin2 (left) and ng = 100 gal/amin2 (right).
4. Multiscale Entropy Restoration 4.1. The Multiscale Entropy The Undecimated Isotropic Wavelet Transform (UIWT) decomposes an n × n image I as a superposition of the form I(k, l) = cJ k,l +
where cJ is a coarse or smooth version of the original image I and wj represents the details of I at scale 2−j (see Starck et al.(Starck et al., 1998; Starck and Murtagh, 2002) for details). Thus, the algorithm outputs J + 1 sub-band arrays of size n × n. We will use an indexing convention such that j = 1 corresponds to the finest scale (high frequencies). The Multiscale Entropy concept (Pantin and Starck, 1996) consists in replacing the standard MEM prior (i.e. the Gull and Skilling entropy) by a wavelet based prior. The entropy is now defined as H(I) =
In this approach, the information content of an image is viewed as sum of information at different scales. The function h defines the amount of information relative to a given wavelet coefficient. Several functions have been proposed for h: – LOG-MSE: The Multiscale Entropy function used in (Pantin and Starck, 1996) (we call it LOG-MSE in the following) is defined by: h(wj,k,l ) =
σj | wj,k,l | (15) 2 [wj,k,l − mj − | wj,k,l | log( K σ )] σX m j
where σX is the total noise standard deviation of the data and σj is the noise standard deviation at scale j. Km is a user-supplied parameter. – ENERGY-MSE: The entropy can be defined as the function of the square of the wavelet coefficients (Starck et al., 2001): h(wj,k,l ) =
2 wj,k,l σj2
The same multiscale entropy function was also derived in Maisinger et al. (2004). – NOISE-MSE: In Starck et al. (2001), the entropy is derived using a modeling of the noise contained in the data: Z |wj,k,l | ∂h(x) h(wj,k,l ) = Pn (| wj,k,l | −u)( )x=u du(17) ∂x 0 where Pn (wj,k,l ) is the probability that the coefficient wj,k,l can be due to the noise: Pn (wj,k,l ) = Prob(W >| wj,k,l |). For Gaussian noise, we have: Z +∞ 2 exp(−W 2 /2σj2 )dW Pn (wj,k,l ) = √ 2πσj |wj,k,l | | wj,k,l | = erfc( √ ) (18) 2σj and 1 h(wj,k,l ) = 2 σj
| wj,k,l | −u √ )du 2σj
LOG-MSE presents an indeterminacy when the wavelet coefficient is equal or close to 0 and the model used in equation (15) is somewhat ad hoc. This point was raised in Maisinger et al. (2004). A better choice for the LOG-MSE would be the Herbert and Leaby function (Hebert and Leahy, 1989) (see also the discussion in section 4.4): | wj,k,l | h(wj,k,l ) ∝ log 1 + (20) σj The ENERGY-MSE is quadratic and leads to a strong penalization even for wavelet coefficients with high signalto-noise ratio. Such penalization terms are known to oversmooth the strongest peaks and should not be used for the weak lensing mass reconstruction. The NOISE-MSE is very close to the l1 norm (i.e. absolute value of the wavelet coefficient) when the coefficient value is large, which is known to produce good results for the analysis of piecewise smooth images (Donoho and Elad, 2003). We therefore choose the NOISE-MSE entropy as the most appropriate for the weak lensing reconstruction problem. Fig. 5
Weak Lensing Mass Reconstruction
shows the l1 norm penalization function, the ENERGYMSE and NOISE-MSE. The NOISE-MSE penalization presents a quadratic behavior for small coefficients and a linear one for larger coefficients. More details are given in section 4.4.
4.2. Significant Wavelet Coefficients using the FDR In Pantin and Starck (1996), it has been suggested to not apply the regularization on wavelet coefficients which are clearly detected (i.e. significant wavelet coefficients). The new Multiscale Entropy is: ¯ (j, k, l)h(wj,k,l ) hn (wj,k,l ) = M
¯ (j, k, l) = 1 − M (j, k, l), and M is the multiresowhere M lution support (Murtagh et al., 1995): 1 if wj,k,l is significant M (j, k, l) = (22) 0 if wj,k,l is not significant This describes, in a Boolean way, whether the data contains information at a given scale j and at a given position (k, l). Commonly, wj,k,l is said to be significant if the probability that the wavelet coefficient is due to noise is small, i.e. if P (| W > wj,k,l |) < ǫ, where P is a given noise distribution function. In the case of Gaussian noise, this amount to state that wj,k,l is significant if | wj,k,l |> kσj , where σj is the noise standard deviation at scale j, and k is a constant, generally taken between 3 and 5 (Murtagh et al., 1995). With this definition, the number of false detections depends on both the ǫ value and the image size. An alternative approach to this detection strategy is the False Discovery Rate method (FDR) (Benjamini and Hochberg, 1995). This technique has recently be introduced for astronomical data analysis (Miller et al., 2001; Hopkins et al., 2002). It allows us to control the average fraction of false detections made over the total number of detections. It also offers an effective way to select an adaptive threshold. The FDR is given by the ratio : F DR =
is no bigger than α: E(F DR) ≤
Ti .α ≤ α V
The unknown factor TVi is the proportion of truly inactive pixels. A complete description of the FDR method can be found in Miller et al. (2001). In Hopkins et al. (2002), it has been shown that the FDR outperforms standard methods for sources detection. Here, we use the FDR method; at each resolution level j of the decomposition. We derive a detection threshold Tj (from a αj value). We have chosen to take a different α value per scale. To fix a α value per scale, we used The Receiver Operating Characteristic (ROC) (Genovese and Eddy, 1997) curves in order to quantify the quality of the detection at a given scale for different α values. We found that the αj value must increase with scale following the relation: αj = α0 ∗2j for the spatial observations and αj = α0 ∗(1.7)j for the ground observations where α0 = 0.0125. We then consider a wavelet coefficient. wj,k,l as significant if its absolute value is larger than Tj .
4.3. Multiscale Entropy Restoration Assuming Gaussian noise, the Multiscale Entropy restoration method lead to the minimization of the functional, J X X −κ ˜ k2 hn ((W κ ˜ )j,k,l ) + β 2σn2 j=1
J(˜ κ) =
where σn the noise standard deviation in κ ˜ l , J the number of scales, β is the regularization parameter and W is the Wavelet Transform operator. The β parameter is calculated automatically under the constraint that the residual should have a standard deviation equal to the noise standard deviation. Full details of the minimization algorithm can be found in Starck et al. (2001), as well as the way to determine automatically the regularization parameter β.
where Via are the number of pixels truly inactive declared active and Da are the number of pixels declared active. This procedure controlling the FDR specifies a rate α between 0 and 1 and ensures that, on average, the FDR
4.4. Related Work 4.4.1. The Generalized Wavelet Regularization Using a prior such that a pixel value is a function of its neighborhood (see Molina et al. (2001) for more details
on the Markov Random Field model), the Bayesian solution consists in adding the following penalization on the solution: XX C(˜ κ) = β φ(˜ κ(x, y) − κ ˜ (x, y + 1))2 + x
+ φ(˜ κ(x, y) − κ ˜(x + 1, y))2
The function φ, called potentialPfunction, is an edge preP serving function. The term β x y φ(k ∇I k (x, y)) can also be interpreted as the Gibbs energy of a Markov Random Field. Generally, functions φ are chosen with a quadratic part which ensures a good smoothing of small gradients (Green, 1990), and a linear behavior which cancels the penalization of large gradients (Bouman and Sauer, 1993): ′
Such functions are often called L2 -L1 functions. Examples of φ functions: 1. φq (x) = x2 : quadratic function. 2. φT V (x) =| x |: Total√Variation. Hyper-Surface 3. φ2 (x) = 2 1 + x2 − 2: (Charbonnier et al., 1997). 4. φ3 (x) = x2 /(1 + x2 ) (Geman and McClure, 1985). 2 5. φ4 (x) = 1 − e−x (Perona and Malik, 1990). 6. φ5 (x) = log(1 + x2 ) (Hebert and Leahy, 1989). Figure 4 shows different φ functions. It has been shown that this concept can be generalized in the wavelet domain, leading to a multiscale wavelet penalization term (Jalobeanu, 2001): X Cw (˜ κ) = β φ(k (W κ ˜ )j,k,l kp ) (26) j,k,l
When φ(x) = x and p = 1, it corresponds to the l1 norm of the wavelet coefficients. In this framework, the multiscale entropy deconvolution method is only one special case of the wavelet constraint deconvolution method. Figure 5 shows the multiscale entropy penalization function. The dashed line corresponds to a l1 penalization (i.e. φ(w) =| w |), the dotted line to a l2 penaliza2 tion φ(w) = w2 , and the continuous line to the multiscale
entropy function. We can immediately see that the multiscale entropy function presents a quadratic behavior for small values, and is closer to the l1 penalization function for large values. Penalization function with a l2 -l1 behavior are known to be a good choice for image restoration.
4.4.2. Multiscale MEM and ICF The multichannel ICF-MEM method (Weir, 1991, 1992) consists in assuming that the visible-space image O is formed by a weighted sum of the visible-space image PNc channels Oj , O = j=1 pj Oj where Nc is the number of channels and Oj is the result of the convolution between a hidden image hj with a low-pass filter (ICF) Cj , called ICF (Intrinsic Correlation Function) (i.e. Oj = Cj ∗hj ). In practice, the ICF is a Gaussian. The MEM-ICF constraint is: Nc X | hj | (27) | hj | −mj − | hj | log CICF = mj j=1 In Maisinger et al. (2004), it was argued that the multiscale entropy is merely a special case of the intrinsic correlation function approach, where we replace the ICF kernel by a wavelet function. From the strict mathematical point of view, this is right, but this vision minimizes completely the improvement related to the wavelets. All the concepts of sparse representation (which is the key of the wavelet success in many applications), fast decomposition and reconstruction, zero mean coefficients (which allows us to get wavelet coefficients which are independent of the background and to derive a robust noise modeling) do not exist in the ICF-MEM approach. Furthermore, ICF-MEM approach requires to estimate accurately the background, which may be sometimes a very difficult task, and it has be shown (Bontekoe et al., 1994) that the solution depends strongly on this estimation. On the contrary, Multiscale MEM needs only an estimation of the noise standard deviation, which is easy to determine. For all these reasons, we prefer to keep our vision of the multiscale entropy method as a specific case of the
Weak Lensing Mass Reconstruction
Figure 4. Examples of potential function φ. generalized wavelet regularization techniques rather than as an extension of the ICF approach.
5. Results 5.1. Comparison of methods We have used a simulated data set obtained using a standard Λ-CDM cosmological model. A part of the κ mass map and the shear maps is shown in Fig. 1. The field size is 2 × 2 square degrees, sampled with 1024 ∗ 1024 pixels. Noisy shear maps, corresponding to both spatial (i.e. ng = 100 gals amin−2 ) and ground-based observations (i.e. ng = 20 gals amin−2 ), are created using equation 7. Then we have reconstructed the two noisy mass maps from equation 9 and applied the following methods: 1. Gaussian filtering with a standard deviation equal to σG = 1 amin. 2. Gaussian filtering with a standard deviation equal to σG = 2.5 amin. 3. Wiener filtering. 4. Maximum Entropy Method (MEM) using the LensEnt2 package. As this code has not been designed for manipulating large images, we had to restrict the restoration by this method to a field size of0.5 × 0.5 square degree, sampled with 256*256 pixels. Since the LensEnt2 maps are positivity constrained, as recommended by the author of the LensEnt2 package, we have recovered a physical mass by transforming the outputs such that the minimum convergence in the central quarter of the reconstruction is zero. To optimize the ICF, we have maximized the Bayesian evidence value as a function of ICF width, and found that maximum evidence is around 210 arcsec for ground observations and around 180 arcsec for space observations. 5. Multiscale Entropy method. The evaluation is done by i) visual inspection of the images, ii) calculating the standard deviation between the original κ mass map and the reconstructed map (i.e. E = ST D(κ−˜ κ) ST D(κ) ), iii) calculating the standard deviation for each of their wavelet scales (i.e. Ej =
ST D((Wκ)j −(W κ ˜ )j ) ) ST D((Wκ)j )
iv) calculation the power spectrum of the error E (for MEM and multiscale entropy methods). The values σG = 1 amin and σG = 2.5 amin have been chosen to optimize the Gaussian filtering for ng = 100 gals arcmin−2 and ng = 20 gals arcmin−2 , respectively. Table 1 gives the standard deviation of the error for the four reconstructed mass maps. It shows that i) the Wiener filtering is better than the Gaussian filtering and the MEM-LensEnt2 method and ii) the Multiscale Entropy outperforms the three other methods. Fig. 6 shows the error versus the scale (each wavelet scale) for both simulations using the Gaussian filtering (continuous line), the Wiener filtering (dotted line), the MEM-LensEnt2 filtering (dashed line) and the Multiscale Entropy filtering (dotted-dashed line). The wavelet scales 1 to 6 correspond to scales of 0.12, 0.23, 0.47, 0.94, 1.87, 3.75 amin respectively. We can see that the Multiscale Entropy method produces better results for all scales. Fig. 7 shows the log power spectrum of the error. It is very consistent with the previous one. Indeed, the MEM error becomes very important toward the smallest frequencies (largest wavelet scales). The same experiment has been done with a smallest ICF (ICF=120 for the spatial simulation), but the result is worse, which is not surprising since the ICF value was chosen to get the best results. Fig. 8 shows from top to bottom the reconstructed maps for the Gaussian, the Wiener and Multiscale Entropy filtering. Fig. 8 left corresponds to ground-based observations (i.e. ng = 20) and Fig. 8 right corresponds to spatial observations (i.e. ng = 100). Fig. 9 shows the denoising results on a portion of the previous image. Fig. 9 shows the original noise free simulated image of the 0.5 × 0.5 square degrees field (upper left), the Multiscale Entropy Filtering for the spatial simulated observations (nG = 100) (upper right), the MEMLensEnt2 restoration for the ground based observations (bottom left) and the spatial observations (bottom right). The computation time for the 1024*1024 pixels map is 4 minutes for the Multiscale Entropy method, 26 seconds for the Wiener filtering and 4 seconds for the Gaussian smoothing. The computation time for the 256*256 pixels map is around 60 minutes (it depends on the convergence of the result) using the MEM-LensEnt2 package.
Table 1. Standard deviation of the reconstruction error with five different methods.
Figure 6. Standard deviation versus scale for the ground-based simulation (left) and the space-based simulation (right). Figure 7. Log Power Spectrum of the Error by Multiscale Entropy Filter and MEM for the ground-based simulation (left) and the space-based simulation (right).
5.2. Robustness to missing data
5.3. Cluster detection
During the observations, various problem can cause a loss of data in the image. For example, it can be due to a defect of the camera CCD, generating a dark line or a dark row in the image, or to the presence of a bright star in the field of view which forces us to remove part of the image. In order to study this problem, we mask two rectangular areas, setting all pixel values to 0, in the shear maps γ1 and γ2 . By inverse filtering, we have derived the noisy mass map κl in which we can also visualize the lack of data (Fig. 10 upper left). Then we have applied the three methods, Gaussian filtering, Wiener filtering and Multiscale Entropy, to the noisy mass map and the results can be seen respectively in Fig. 10 upper right, Fig. 10 bottom left and bottom right. We can see that all three methods are robust to the missing data. Note however that, for the Wiener filtering, we have assumed perfect knowledge of the power spectrum of κ, while, in practice, its estimation is made more complicated by the complex field geometry.
Another important aspect of the weak shear mass reconstruction is the possibility to detect clusters and to build a catalog. Here, using the FDR in the wavelet space, we detect as significant a set of wavelet coefficients. We built an isophote map, where each isophote level corresponds to the detection level in a given scale. This isophote is overplotted on the true mass map, which allows us to visually check the false detections and the missed detections. A cluster surrounded by two isophotes means that it has been detected at two scales. Fig. 12 left shows the isophote map when we use the regular kσ thresholding and Fig. 12 left right shows the isophote map when we use the FDR method. We see that the FDR is more sensitive than the kσ method for the detection, without being contaminated by a large number of false detections. Fig. 13 shows a zoom of these two maps. Figure 14 shows a comparison between the Gaussian filtering, the Wiener filtering and FDR-Wavelet method for the detection of clusters. In the Gaussian and Wiener maps, the isophotes corresponds to a kσ detection level where k = 3, 4, 5. It shows clearly how the FDR-Wavelet method outperforms the other methods.
Fig. 11 shows the error versus the scale for both simulations using the Gaussian filtering (continuous line), the Wiener filtering (dotted line) and Multiscale Entropy. We can see that the Multiscale Entropy still produces better results at all scales. Bayesian methods such MEM could also take into account properly missing data, however not in a straightforward way as when using wavelets.
5.4. E/B Decomposition As explained in §2.3, a simple diagnostic test for a wide range of systematic effects is to search for the presence of B-mode in the lensing maps. In order to test it, we have simulated mass maps with a B-mode.
Weak Lensing Mass Reconstruction
Figure 8. Restoration of the 2 × 2 square degrees ground-based observation (left) and spatial observation (right). From top to bottom, Gaussian filtering, Wiener filtering and Multiscale Entropy filtering.
Figure 9. In a region of 0.5 × 0.5 square degrees, a sixteenth of the original field : Upper left, simulated mass map, upper right, Multiscale entropy filtering for ng = 100 gal/arcmin2. Bottom left, MEM filtering for ng = 20 gals/amin−2 (ICF width = 210) and bottom right for ng = 100 gals/amin−2 (ICF width = 180).
Figure 16. Noisy simulated mass map Fig. 15 left shows a simulated mass map with a lensing E-mode signal (left) and an arbitrary B-mode signal (left). As usual, we have added a realistic space-based Gaussian noise to the shear of this simulation. Fig. 16 shows the noisy mass map resulting. Using the Multiscale Entropy filtering, we have then reconstructed the two components of the mass map (see §2.3): E-mode in Fig. 17 left and B-mode in Fig. 17 right. We see clearly that the wavelet separation of the E and B modes is very good. Indeed, the two main features in the B-mode have well been recovered, without interfering with the reconstruction of the E-mode.
6. Conclusion We have presented in this paper a new way to reconstruct weak lensing mass maps. We have modified the Multiscale Entropy method in order to take into account the FDR. We have shown that this new method outperforms several standard techniques currently used for the weak shear mass reconstruction. The visual aspect as well as objective criteria, such the rms of the error or the rms per scale of the error, clearly show the advantages of the proposed approach. Experiments have demonstrated that it is also robust to missing data. We have also shown that a E/B mode separation can also be performed using this method, thus providing a useful test for the spatial distribution of residual systematics. Our method allows us also
to build a catalog of clusters and the use of FDR leads to a clear improvement in sensitivity, compared to what has been done previously with wavelets.
Software The software related to this paper, MR/Lens, and its full documentation are available from the following web page:
Acknowledgements. We wish to thank Savita Mathur for her initial work on the subject, Richard Massey and Bedros Afeyan for useful discussions. We would like to particularly acknowledge Phil Marshall for his help relative to the use of the LensEnt2 package.
References Bartelmann, M.: 1995, A&A 303, 643 Bartelmann, M., Narayan, R., Seitz, S., and Schneider, P.: 1996, ApJ 464, L115+ Bartelmann, M. and Schneider, P.: 1999, ArXiv Astrophysics e-prints Benjamini, Y. and Hochberg, Y.: 1995, J. R. Stat. Soc. B 57, 289 Bontekoe, T., Koper, E., and Kester, D.: 1994, A&A 284, 1037
Figure 11. Standard deviation versus scale for the ground-based simulation(left) and the spatial simulation (right) with missing data Figure 12. Isophotes of the detected wavelet coefficients for the space-based simulation overplotted on the original mass map: left, using the kσ standard approach (with k = 3, 4, 5) and right using the FDR method. Figure 13. Zoom of the previous maps Bouman, C. A. and Sauer, K.: 1993, IEEE Transactions on Image Processing 2(3), 296 Bridle, S. L., Hobson, M. P., Lasenby, A. N., and Saunders, R.: 1998, MNRAS 299, 895 Charbonnier, P., Blanc-Feraud, L., Aubert, G., and Barlaud, M.: 1997, IEEE Transactions on Image Processing 6(2), 298 Donoho, D. L. and Elad, M.: 2003, the Proc. Nat. Aca. Sci. 100, 2197 Geman, S. and McClure, D.: 1985, in A. S. Assoc. (ed.), Proc. Statist. Comput. Sect., Washington DC Genovese, C. R.and Noll, N. C. and Eddy, W. F.: 1997, Magn. Res. Med. 38, 497 Green, P. J.: 1990, IEEE Transactions on Medical Imaging 9(1), 84 Gull, S. and Skilling, J.: 1991, MEMSYS5 Quantified Maximum Entropy User’s Manual, Royston,England Hebert, T. and Leahy, R.: 1989, IEEE Transactions on Medical Imaging 8(2), 194 Hopkins, A. M., Miller, C. J., Connolly, A. J., Genovese, C., Nichol, R. C., and Wasserman, L.: 2002, AJ 123, 1086 Jalobeanu, A.: 2001, Ph.D. thesis, Universit´e de Nice Sophia Antipolis Kaiser, N.: 1995, ApJ 439, L1 Kaiser, N. and Squires, G.: 1993, ApJ 404, 441 Maisinger, K., Hobson, M. P., and Lasenby, A. N.: 2004, MNRAS 347, 339 Marshall, P. J., Hobson, M. P., Gull, S. F., and Bridle, S. L.: 2002, MNRAS 335, 1037 Mellier, Y.: 1999, ARA&A 37, 127 Mellier, Y.: 2002, Space Science Reviews 100, 73 Miller, C. J., Genovese, C., Nichol, R. C., Wasserman, L., Connolly, A., Reichart, D., Hopkins, A., Schneider, J., and Moore, A.: 2001, AJ 122, 3492 Molina, R., N´ un ˜ ez, J., Cortijo, F., and Mateos, J.: 2001, IEEE Signal Processing Magazine 18, 11
Murtagh, F., Starck, J.-L., and Bijaoui, A.: 1995, A&AS 112, 179 Narayan, R. and Nityananda, R.: 1986, ARA&A 24, 127 Pantin, E. and Starck, J.-L.: 1996, A&AS 315, 575 Perona, P. and Malik, J.: 1990, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(7), 629 Refregier, A.: 2003, ARA&A 41, 645 Schneider, P. and Seitz, C.: 1995, A&A 294, 411 Seitz, S., Schneider, P., and Bartelmann, M.: 1998, A&A 337, 325 Squires, G. and Kaiser, N.: 1996, ApJ 473, 65 Starck, J.-L. and Murtagh, F.: 2002, Astronomical Image and Data Analysis, Springer-Verlag Starck, J.-L., Murtagh, F., and Bijaoui, A.: 1998, Image Processing and Data Analysis: The Multiscale Approach, Cambridge University Press Starck, J.-L., Murtagh, F., Querre, P., and Bonnarel, F.: 2001, A&A 368, 730 Starck, J.-L., Pantin, E., and Murtagh, F.: 2002, Publications of the Astronomical Society of the Pacific 114, 1051 Vale, C. and White, M.: 2003, ApJ 592, 699 Van Waerbeke, L., Mellier, Y., Radovich, M., Bertin, E., Dantel-Fort, M., McCracken, H. J., Le F`evre, O., Foucaud, S., Cuillandre, J.-C., Erben, T., Jain, B., Schneider, P., Bernardeau, F., and Fort, B.: 2001, A&A 374, 757 Weir, N.: 1991, in 3rd ESO/ST-ECF Data Analysis Workshop Weir, N.: 1992, in D. Worral, C. Biemesderfer, and J. Barnes (eds.), Astronomical Data Analysis Software and System 1, pp 186–190, Astronomical Society of the Pacific
Weak Lensing Mass Reconstruction
Figure 14. The isophotes represent the detected clusters using the Gaussian filtering (upper left), the Wiener filtering (upper right) and the wavelet-FDR method. Figure 15. left mass map (E-mode), right mass map (B-mode) Figure 17. left filtered noisy mass map (E-mode), right filtered noisy mass map (B-mode)