Improving Weak Lensing Mass Map Reconstructions using Gaussian and Sparsity Priors: Application to DES SV
N. Jeffrey F. B. Abdalla O. Lahav F. Lanusse J. -L. Starck A. Leonard D. Kirk C. Chang E. Baxter T. Kacprzak S. Seitz V. Vikram L. Whiteway T. M. C. Abbott S. Allam S. Avila E. Bertin D. Brooks A. Carnero Rosell M. Carrasco Kind J. Carretero F. J. Castander M. Crocce C. E. Cunha C. B. D'Andrea L. N. da Costa C. Davis J. De Vicente S. Desai P. Doel T. F. Eifler A. E. Evrard B. Flaugher P. Fosalba J. Frieman J. Garcia-Bellido D. W. Gerdes D. Gruen R. A. Gruendl J. Gschwend G. Gutierrez W. G. Hartley K. Honscheid B. Hoyle D. J. James M. Jarvis K. Kuehn M. Lima H. Lin M. March P. Melchior F. Menanteau R. Miquel A. A. Plazas K. Reil A. Roodman E. Sanchez V. Scarpine M. Schubnell I. Sevilla-Noarbe M. Smith M. Soares-Santos F. Sobreira E. Suchyta M. E. C. Swanson G. Tarle D. Thomas A. R. Walker
astro-ph.CO


firstpage–lastpage
Out-of-time-ordered measurements as a probe of quantum dynamics
    Pavan Hosur
    December 30, 2023
===============================================================





Mapping the underlying density field, including non-visible dark matter, using weak gravitational lensing measurements is now a standard tool in cosmology. Due to its importance to the science results of current and upcoming surveys, the quality of the convergence reconstruction methods should be well understood. We compare three methods: Kaiser-Squires (KS), Wiener filter, and Glimpse. KS is a direct inversion, not accounting for survey masks or noise. The Wiener filter is well-motivated for Gaussian density fields in a Bayesian framework. Glimpse uses sparsity, aiming to reconstruct non-linearities in the density field. We compare these methods with several tests using public Dark Energy Survey (DES) Science Verification (SV) data and realistic DES simulations. The Wiener filter and Glimpse offer substantial improvements over smoothed KS with a range of metrics. Both the Wiener filter and Glimpse convergence reconstructions show a 12 per cent improvement in Pearson correlation with the underlying truth from simulations. To compare the mapping methods' abilities to find mass peaks, we measure the difference between peak counts from simulated ΛCDM shear catalogues and catalogues with no mass fluctuations (a standard data vector when inferring cosmology from peak statistics); the maximum signal-to-noise of these peak statistics is increased by a factor of 3.5 for the Wiener filter and 9 for Glimpse. With simulations we measure the reconstruction of the harmonic phases; the phase residuals' concentration is improved 17 per cent by Glimpse and 18 per cent by the Wiener filter. The correlation between reconstructions from data and foreground redMaPPer clusters is increased 18 per cent by the Wiener filter and 32 per cent by Glimpse. 





gravitational lensing: weak  – large-scale structure of Universe– methods: statistical 









§ INTRODUCTION
 

Mass map reconstruction from weak gravitational lensing recovers the underlying matter distribution in the universe from measurements of galaxy shapes. Images of distant galaxies are deformed by the inhomogeneous matter distribution along the line of sight. Any matter can contribute to the lensing effect, making it a direct probe of non-visible dark matter.

Weak lensing, which takes advantage of the statistical power from many small distortions (that is, gravitational lensing induced “shears”), is now a well established tool in constraining cosmology. The Dark Energy Survey (DES) has used the 2-point correlation function of shear to contribute to excellent constraints on cosmological parameters and models, including the nature of dark energy <cit.>. Shear 2-point correlation functions have been used to constrain cosmology from many other survey datasets (, ). These methods use the shear measurements directly, as the shear can be related to the underlying matter distribution without needing to explicitly reconstruct mass maps.

A zero-mean Gaussian random field can be characterised entirely by its 2-point correlations. The matter density field in the early universe is expected to be highly Gaussian, a property which persists into the late universe for the large scales that were less affected by gravitational collapse. For the smaller scales at late times, non-linear gravitational collapse has led to a highly non-Gaussian density field. Much valuable information can be extracted from this non-Gaussianity, although this requires additional methods beyond 2-point statistics.

Popular proposed methods to extract this information include N-point statistics and higher order moments <cit.>, peak statistics (, , , , ), and Minkowksi functionals (, ). It is often either essential or convenient to apply these methods to the density field directly (rather than in the space of the shear measurements), thereby necessitating a reliable mass map reconstruction. 

Peak statistics are particularly promising, as peaks in the density field probe the non-Gaussian structure directly. Peaks can be identified from aperture mass maps, which are derived by convolving the shear data with a kernel, or from the reconstructed density field. The first approach has the advantage of having local noise, while the second is “closer” to the underlying density field and often has faster algorithms. Both methods often require simulations to provide a link between the theory and data, with the exception of proposed semi-analytic models (, ).

In addition to using mass maps for higher order statistics to constrain cosmological parameters and models, the mass maps can themselves be intrinsically useful. <cit.>, using the original DES Science Verification (SV) mass map, show evidence that the 1-point distribution of the density field is more consistent with Log-Normal than Gaussian. Combining mass maps with the spatial distributions of stellar mass, galaxies, or galaxy clusters allows the relationship between the visible baryonic matter and invisible dark matter to be studied. Using mass maps to constrain galaxy bias <cit.>, the relation between the distribution of galaxies and matter, can in turn aid cosmological probes other than weak lensing. Maps also enable simple tests for systematic error in the galaxy shape catalogues.

Since the first application of mass mapping methods to wide-field surveys with the Canada-France Hawaii Telescope Lensing Survey (CFHTLenS) data <cit.>, mass maps have been a standard product of large weak lensing surveys. In addition to DES, current surveys reconstructing the density field from weak lensing data include the Kilo-Degree Survey (Giblin et al. in prep.) and the Hyper Supreme-Cam Subaru Strategic Program (HSC-SSP) <cit.>. Mapping dark matter is key to the science goals of the future Euclid Mission <cit.> and the Large Synoptic Survey Telescope <cit.>.

DES is a ground based photometric galaxy survey, observing in the southern sky from the 4m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO) in Chile with five photometric filters covering the optical and near-infrared spectrum using the Dark Energy Camera (, ). The SV data come from an initial run over a fraction of the final sky coverage, but to almost the full exposure time of the final survey. The sky coverage is still large, 139 deg^2, and the nearly full exposure <cit.> gives a galaxy density almost equal to what is expected after the complete 5 years of DES observations.

This paper uses the public DES SV data to compare the quality of mass mapping reconstruction methods. The maps are of the two-dimensional convergence, κ, a weighted projection of the density field in the foreground of the observed background galaxies. Recovering the convergence from the shear data is an ill-posed inverse problem, troubled by survey masks and galaxy “shape noise”. 

This work follows on from that of <cit.> and <cit.> in which the original DES SV mass map was created using the <cit.> method. In this paper we compare three quite different methods: Kaiser-Squires (KS); Wiener filtering <cit.>; and Glimpse (, ), a sparsity-based reconstruction method. The Kaiser-Squires method is a direct inversion from shear to convergence, taking no account of missing data or the effect of noise. The Wiener filter and Glimpse assume different prior knowledge about the underlying convergence to account for the effects of noise and missing data.

In Sec. 2 we describe the theoretical foundation for weak lensing mass mapping and the three different methods used for this work. In Sec. 3 we describe the DES SV shear data, the accompanying simulations, and the redMaPPer galaxy cluster catalogue. Foreground galaxy clusters are expected to trace the true density field, and therefore should be correlated with the convergence reconstruction. The different methods are also applied to realistic data simulations where the true convergence is known. In Sec. 4 we present our results on data and simulation, using various quality metrics for the reconstruction. On simulations these metrics are the Pearson correlation coefficient, the pixel root-mean-square error (RMSE), the variance of the 1-point distribution of pixel values, the phase residuals, and peak statistics. On data we compare the convergence reconstructions to the foreground galaxy clusters. We conclude in Sec. 5.



§ METHODOLOGY



 §.§ Weak gravitational lensing


We can use measurements of the distortion of background galaxy shapes by weak gravitational lensing to learn about the mass distribution in the foreground without making many physical assumptions or relying on phenomenological models. For convenience, here we summarise some of the existing literature relevant for mass mapping from weak lensing (, ). 

The weak lensing formalism follows photon paths along geodesics in a perturbed Friedmann-Robertson-Walker (FRW) metric. The perturbations are sourced by the density field of large scale structure. Throughout we assume that the perturbations are small, and that the measurements are made over a small enough patch of the sky that the sky geometry is Euclidian. Consistent with the Planck CMB results <cit.> and motivated by inflationary theory, we assume that the global geometry of the universe is flat.

The density contrast, δ = (ρ - ρ̅)/ρ̅, of a pressureless fluid is related to the scalar gravitational potential perturbation, Φ, through the Poisson equation,


    ∇^2 Φ = 3 H_0 Ω_m/2 a δ ,


where H_0 is the present value of the Hubble parameter, a is the cosmological scale factor, and ρ and ρ̅ are the local and mean density respectively.

For a flat universe, the lensing potential is given by


    ψ(θ, ω) = 2 ∫_0^ωdω' [ ω - ω'/ωω'] Φ(θ , ω')  ,


where ω is the comoving distance.

The Born approximation assumes that the observed angle to a point, θ, deviates only a small amount from the true angle β, so the change in distance of the photon's path is negligible. We can characterise the effect of lensing on the galaxies using the Jacobian of the transformation, 𝒜_ij = ∂β_i / ∂θ_j, which is decomposed into the functions κ(θ) and γ(θ) = γ_1 + i γ_2, and which is given by 


    𝒜   = ( δ_ij -  ∂^2 ψ (θ)/∂θ_i ∂θ_j )
    
       = ( [ 1 - κ - γ_1        -γ_2;        -γ_2 1 - κ + γ_1 ])  .


Using the definition of the lensing potential and the Poisson equation, the convergence can be expressed as an integral over the density along the line of sight,


    κ(θ, ω ) = 3 H_0^2 Ω_m/2∫_0^ωdω' ω' (ω - ω')/ωδ(θ, ω')/a(ω') .


For a distribution n(ω) of lensed galaxies, the lensing efficiency kernel is defined to be


    p(ω') = ∫^ω_∞_ω^'( ω - ω^'/ω)n(ω) dω ;


this weights the contribution of the foreground density fluctuations to give the convergence weighted over the redshift distribution of source galaxies, 


    κ (θ)    = ∫_0^∞ n(ω) κ(θ, ω ) dω
       = 3 H_0^2 Ω_m/2∫_0^∞dω' p(ω')  ω' δ(θ, ω')/a(ω') .


The shear, γ( θ ), which is assumed to be an observable in the weak lensing limit, is given by


    γ( θ ) = 1/π∫_ℝ^2d^2 θ' 𝒟 ( θ - θ') κ(θ') 
    where 𝒟 ( θ ) = - ( θ_1 - i θ_2 ) ^-2 .


For surveys where the integral is over large angles on the sky, this formulation breaks down, and requires a full treatment in spherical bases. <cit.> show that errors can be introduced at an 𝒪(1 per cent) level for correlations between points at DES SV angular separation depending on the projection. All of the methods used here use the small angle approximation, and should suffer equally.

The real and imaginary parts of the shear γ represent a chosen two dimensional coordinate system. In weak lensing, the observed ellipticity[Using the <cit.> equation 4.10 ellipticity definition for ϵ.] of a galaxy ϵ_obs is related to the reduced shear g plus the intrinsic ellipticity of the source galaxy ϵ_s through


    ϵ_obs   ≈ g + ϵ_s 
         where  g = γ/1 - κ .


The reduced shear is approximately the true shear, g ≈γ, in the weak lensing limit. This allows a standard definition of observed shear, γ_obs = ϵ_obs, where the measurements are degraded by “shape noise”, caused by the ϵ_s values of the observed galaxies:


    γ_obs≈γ +   ϵ_s  .


The shape noise for a given galaxy is modelled as a randomly-drawn Gaussian variate, ϵ_s ∼ G(0,σ_ϵ), where σ_ϵ is estimated from data. The distribution of the ellipticity from the SV data in figure <ref> is not an exact Gaussian, as the true distribution is the result of galaxy astrophysics, though a Gaussian still has properties that make it a good approximation. The Gaussian would be the maximum entropic, least informative, distribution for known mean and variance, and, by the central limit theorem, would be the correct distribution in the limit of large numbers of galaxies averaged in pixels.

It is possible to extend the simple Kaiser-Squires method (Sec. <ref>) to use the reduced shear, g, for the mildly non-linear lensing regime when it is no longer appropriate to assume g ≈γ (, , ). This is also done by Glimpse (Sec. <ref>).

In matrix notation, the problem as given by equations <ref> and <ref> can be expressed as a linear model, with a data vector of observed shear measurements 


    γ = 𝐀κ + 𝐧 ,


where 𝐀 is a discretised version of equation <ref> and 𝐧 is a noise vector due to shape noise (equation <ref>). The elements of the data vector can either correspond to the individual shear measurements or to measurements binned into angular pixels (in which case the noise vector would be the average noise in the pixel). 

The convergence need not be reconstructed with the same pixelisation as the shear measurements, giving κ and γ vectors of different length. Missing data due to survey masks would correspond to a shorter γ vector; here one may wish to fill in the convergence in the masked region — this is known as inpainting. Different sized  κ and γ vectors result in a  non-square 𝐀 matrix, potentially causing inversion problems.



 §.§ Kaiser-Squires reconstruction
 


  §.§.§ Theory


The convergence-to-shear relationship, equation <ref>, is a convolution in the two dimensional angular plane. The two-dimensional Fourier transforms of the shear and convergence, defined for κ as


    κ̃ ( k ) = ∫_ℝ^2d^2 θκ(θ) exp ( i θ·k )  ,


are related through an elementwise product via the convolution theorem


    γ̃ ( k ) = π^-1𝒟̃ ( k ) κ̃ ( k )  ,


where the Fourier transform of the kernel is given by


    𝒟̃ ( k ) = π(k_1^2 - k_2^2 + 2 i k_1 k_2) /|k|^2 ;


here k_1 and k_2 are the components of k. Using 𝒟̃𝒟̃^* = π^2, equation <ref> can be rewritten:


    κ̃ ( k ) = π^-1γ̃ ( k )  𝒟̃^* ( k )     for  k≠0 .


The inverse Fourier transform then returns the convergence reconstruction in configuration space <cit.>. 

The real and imaginary parts of the reconstruction are the E- and B-modes respectively, where κ_recon = κ_E + iκ_B. In standard cosmology (equation <ref>), the convergence sourced by a real density field should be a pure E-mode. Errors, noise or other systematic effects can lead to B-mode contributions to the reconstruction. 



  §.§.§ Implementation


In the matrix formulation of equation <ref>, this deconvolution corresponds to multiplying the Fourier space shear field with the inverse of 𝐀 in Fourier space. For a case with no shape noise, that is


    γ̃ = 𝐀̃κ̃ ,


the Kaiser-Squires method is identical to using the inverse matrix


    [ 𝐀̃^-1]_ij   =  k_1,i^2 - k_2,i^2 - 2 i  k_1,i k_2,i/k_1,i^2 + k_2,i^2δ_ij
       = [ 𝐀̃^†]_ij ,


where the Kronecker delta function, δ_ij, relates the element-wise multiplication in Fourier space to a diagonal matrix operator, and ^† is the conjugate transpose.

For the Kaiser-Squires inversion in configuration space, the 𝐀 and 𝐀^† matrices are not diagonal, and therefore are slower to compute. The discretisation of the underlying smooth shear field into finite configuration space makes the property 𝐀𝐀^†  = 𝐈 inexact. As a result of these factors, we choose to implement the Kaiser-Squires reconstruction in Fourier space.

The shear due to lensing is much smaller than the shape noise, and not all places on the sky contain usable galaxies. Both the shape noise and the random sampling of background galaxies propagate error through this noisy reconstruction. Binning the shear measurements into larger pixels can reduce the shape noise per pixel and ensure that there are no empty pixels, but this comes at a loss of the small scale information and cannot deal with masks or the edges of the survey.

A smoothing filter is applied to the Kaiser-Squires reconstruction to reduce the noise. This will similarly lose any small scale structure, and especially suppress peaks in the convergence. In this work, matching <cit.>, we smooth the Kaiser-Squires maps with a Gaussian kernel. The standard deviation scale, σ_smooth, of this Gaussian kernel is free to be chosen, where σ_smooth = 0 corresponds to standard, unsmoothed Kaiser-Squires.



 §.§ Wiener Filter



  §.§.§ Theory
 

The Wiener filter is the linear minimum-variance solution to linear problems of the type in equation <ref>, where the noise is uncorrelated. The Wiener filter reconstruction (, ) is given by 


    κ_W    = 𝐖γ
    𝐖   = 𝐒_κ𝐀^†[ 𝐀𝐒_κ𝐀^† + 𝐍]^-1 .


Here 𝐒_κ and 𝐍 are the signal and noise covariance matrices respectively, which are ⟨κκ^†⟩ and ⟨𝐧𝐧^†⟩ for this problem. 

This filter is the linear minimum-variance solution, as  𝐖 is a linear operator that minimises the variance


    ⟨ (𝐖 γ - κ)^† (𝐖 γ - κ) ⟩ .


If the chosen prior on κ does not constrain the reconstruction, so that 𝐒_κ^-1→ 0 <cit.>, or if the data are noise-free, 𝐍 = 0, then the linear minimum variance filter becomes the Kaiser-Squires reconstruction. Setting 𝐒_κ^-1→ 0 is equivalent to removing the signal prior in the following Bayesian framework.

From a different starting point, for the Wiener posterior we begin by assuming a Gaussian likelihood <cit.>


    Pr( γ | κ ) = 1/√((det 2 π𝐍))exp[ - 1/2 ( γ - 𝐀κ )^† 𝐍^-1 ( γ - 𝐀κ )  ]  ,


where it is assumed that 𝐍 is known and the noise is both uncorrelated and Gaussian, as assumed in equation <ref>. Intrinsic alignments of clustered galaxies will violate this uncorrelation condition.

The prior on the convergence is that of a Gaussian random field, which is applicable for the density field on large scales at late times, 


    Pr( κ | 𝐒_κ ) = 1/√((det 2 π𝐒_κ))exp[ - 1/2κ^† 𝐒_κ^-1κ]  .


Using Bayes' theorem and the fact that Pr(γ | 𝐒_κ, κ ) = Pr(γ | κ ), the full posterior is given by


    Pr( κ | 𝐒_κ, γ )    =   Pr( γ | κ ) Pr( κ | 𝐒_κ ) /Pr( γ )
       ∝1/√((det 2 π𝐒_κ))1/√((det 2 π𝐍))×
           exp[ - 1/2κ^†𝐒_κ^-1κ  - 1/2 ( γ - 𝐀κ )^†𝐍^-1 ( γ - 𝐀κ ) ] 
       ∝exp[  - 1/2 ( κ - 𝐖γ )^† (  𝐒_κ^-1 + 𝐀𝐍^-1𝐀^† ) (κ - 𝐖γ) ]  ,


where 𝐖 is the Wiener filter, so the maximum a posteriori (MAP) solution is that of the Wiener reconstruction.

The choice of Gaussian prior is physically motivated for the large, linear scales of the density field (see section <ref>); alternative prior distributions can be used to give different κ posterior distributions that can be maximised or from which samples can be drawn <cit.>. Recent work by <cit.> proposes the use of a Log-Normal prior distribution. This appears to fit the κ distribution from simulations (figure <ref>) and data <cit.> better than Gaussian, but, unlike the Wiener filter, lacks an analytic MAP solution.  

If the aim of the reconstruction is to infer cosmology from the non-Gaussian component of the density field, the Wiener filter may not be the ideal method for mass map recovery. The small scale modes with less power are often suppressed, losing the peak structure. Qualitatively it can be thought of as either the Gaussian prior being inappropriate or as the linear filter being insufficient.



  §.§.§ Implementation


Using the exact Fourier space property 𝐀̃^-1 = 𝐀̃^† we rewrite equation <ref> as


    κ̃_W    = 𝐀̃^-1𝐒̃_γ[ 𝐒̃_γ + 𝐍̃]^-1γ̃
       = 𝐀̃^-1γ̃_W
       =𝐀̃^†γ̃_W ,


where we have used 𝐀̃𝐒̃_κ𝐀̃^† =  ⟨𝐀̃κ̃κ̃^†𝐀̃^†⟩  = ⟨γ̃γ̃^†⟩ = 𝐒̃_γ. This shows that applying the Wiener filter to the shear to recover γ_W and then applying the Kaiser-Squires inversion in Fourier space is equivalent to directly calculating the Wiener filter of the convergence.

In configuration space, the noise covariance matrix is given by


    [ 𝐍]_ij = 2 σ^2_ϵ/p_iδ_ij ,


where p_i is the galaxy count per pixel. Empty pixels in the masked region have infinite variance, absorbing the mask into a special case of the Wiener filter denoising.

The signal properties for a Gaussian random field are constrained entirely by the mean and the signal covariance matrix, which in harmonic space is identical to the power spectrum. The cosmological principle implies that the angular distribution of a field on the sky is statistically isotropic, so the angular power spectrum, C_ℓ, can contain all the 2-point statistical information. The angular power spectrum of the physical shear E-mode shear signal is defined as


    C_ℓ, E   = 1/2 ℓ + 1∑_m = -ℓ^+ℓ⟨ |a_ℓ m, E|^2⟩
       = C_ℓ, κ ,


where a_ℓ m are the spherical harmonic coefficients and the brackets ⟨ ⟩ average over realisations of the signal. The second equality assumes the flat sky approximation for high[We omit a prefactor which goes as 1 - 𝒪(ℓ^-2) for high ℓ.]ℓ.

We generate a theoretical power spectrum using the Limber approximation with the Cosmosis package <cit.> with our prior fiducial cosmological parameters: Ω_m = 0.286, Ω_Λ = 0.714, Ω_b = 0.047, h = 0.7, σ_8 = 0.82, n_s = 0.962  and  w = -1. We use a background galaxy distribution defined from equation <ref>, and shown in figure <ref>.

It is commonly asked whether it is reasonable to assume cosmological parameters in the map reconstruction, if the maps are then used to infer cosmological parameters. Though we assume a specific set of cosmological parameters, it would still be possible to use the maps for cosmological parameter estimation, from peak statistics for example, if the same prior is used on the simulations and the data identically. If simulations are not used, the power spectrum can be jointly inferred from the data <cit.> using Gibbs sampling.

In order to generate the power spectrum in flat Fourier space, rather than on the curved sky, we again use a flat sky approximation


    k^2_θ P(k_θ) =  ( 𝒩/2 π)^2 ℓ ( ℓ + 1) C_ℓ ,


adapted from <cit.>, where 𝒩 is the total number of pixels in the map, k_θ is the magnitude of the projected Fourier mode, and where we have defined our projected angular power spectrum as


    P(k_θ) δ( 𝐤_θ - 𝐤^'_θ)= ⟨ γ̃ (𝐤_θ )  γ̃^†  (𝐤_θ)  ⟩ .


The largest scale mode is ℓ = 20.51, which corresponds to an angular separation of 17.55  deg.

Though the signal covariance matrix is diagonal in harmonic space (equation <ref>), and the independent noise has covariance which is diagonal in configuration space (equation <ref>), there is no natural basis in which both are sparse. Inversion of dense matrices to evaluate the Wiener filter is bypassed using the algorithm presented in <cit.>, where an additional messenger field is used to pass information between harmonic and configuration space, iteratively converging to the Wiener filter solution. 

These messenger field methods were extended by <cit.> to draw Markov chain Monte Carlo (MCMC) samples from the whole Wiener posterior (equation <ref>). The first application of messenger field methods to weak lensing data was by  (, ), who drew samples from the Wiener posterior and generated Wiener filtered shear (not convergence) maps from CFHTLenS data. By comparison, in this work we do not sample from the Wiener posterior; instead, we use the original messenger field algorithm of <cit.> to calculate the Wiener filter reconstruction of the convergence map from DES SV shear data and simulations. 



 §.§ Sparsity reconstruction
 


  §.§.§ Theory

Consider the coefficients α of the decomposition of a signal 𝐱 in a representation space (or “dictionary”) Φ, so that 𝐱 = Φα. Example dictionaries include the Fourier transform or wavelet transforms. Assuming a sparse prior on the signal 𝐱 in the dictionary Φ means that its representation α is expected to be sparse, that is, with most of the coefficients equal to 0 <cit.>. A simple example is a cosine function signal and a Fourier transformation; in this sparse basis only two coefficients have a non-zero value (corresponding to the frequency of the cosine function).

Formally most signals cannot strictly be made sparse, and are merely compressible with a choice of an appropriate transformation, such as a wavelet transform (, ). For a compressible signal the magnitude-ordered sparse coefficients, α_i, are expected to have exponential decay and therefore to have a Laplace distribution <cit.>. 

Consider a generic linear inverse problem of the form 𝐲 = 𝐀𝐱 + 𝐧. A robust estimate of the signal 𝐱 can be recovered by solving the (“LASSO”) optimisation problem


    α   || 𝐲 - 𝐀Φα ||^2_2 + λ || α ||_1   ,


where λ is a Lagrangian multiplier <cit.>. Here the first term corresponds to a χ^2 minimisation, ensuring fidelity of the signal reconstruction, while the second is the sparsity-promoting regularisation term.

We can include non-constant noise variance by weighting the first χ^2 according to the variance. If the noise variance is included in the χ^2 term, the λ value can be interpreted as a signal-to-noise level in the transformed (e.g. wavelet) space.

The second term does not use the Euclidan l_2 norm, but instead uses the sparsity-promoting l_1 norm, defined as


    || α ||_1 = ∑_i | α_i |  .


These methods are non-linear, so it can be difficult to derive properties analytically. With realistic simulations of the data and true signal, the value for λ can be chosen to maximise some success metric. This is analogous to selecting a theoretical power spectrum for the Wiener filter, or a smoothing scale for Kaiser-Squires.

Sparse recovery methods are non-linear and are not necessarily formulated in the Bayesian framework of the Wiener filter. The Wiener filter reconstruction is that which maximises the Wiener posterior, which is known analytically provided the noise and signal are Gaussian with known covariance. However, one may make a frequentist estimate of the error of the sparse reconstruction by propagating the noise properties of the data using bootstrapping or Monte Carlo techniques.



  §.§.§ Implementation/Glimpse


The choice of dictionary depends on the structures contained in the signal. Theory of structure formation in the universe predicts the formation of quasi-spherical halos of bound matter. It is standard practice to represent the spatial distribution of matter in halos with spherically symmetric Navarro-Frenk-White <cit.> or Singular Isothermal Sphere profiles. Coefficients of Isotropic Undecimated Wavelets <cit.> in two dimensions are well suited to the observed convergence of a dark matter halo. The wavelet transform used in the Glimpse algorithm is the starlet <cit.>, which can represent positive, isotropic objects. 

The sparsity prior in the starlet basis enforces a physical model that the matter field is a superposition of spherically symmetric dark matter halos. This is not wholly correct, but is an approximation which is true for the non-linear regime in the standard model of structure formation, similarly to how the assumption of Gaussianity holds in the linear regime. On large scales, where the density field is expected to be Gaussian, the Glimpse sparsity prior is less appropriate.

The Glimpse algorithm aims to solve the optimisation problem


    κ̂ = κ   ||𝐍^-1/2[ γ   - 𝐓^†𝐀̂𝐅κ] ||^2_2 
       + λ || ωΦ^†κ ||_1  +   i_Im(κ) = 0 ,


where 𝐅 is the Fourier transform matrix, 𝐓 is the Non-equispaced Discrete Fourier Transform (NDFT) matrix, 𝐀̂ is defined in equation <ref>, ω is a diagonal matrix of weights, and Φ^† is the inverse wavelet transform. The indicator function i_Im(·) = 0 (defined in Appendix <ref>) in the final term imposes realness on the reconstruction (no B-modes). The use of NDFT allows the first term to perform a forward fitted Kaiser-Squires-like step without binning the shear data, allowing the smaller-scales to be retained in the reconstruction. The full algorithm, including the calculation of the weights, is described in Sec. 3.2 in <cit.>.

Though the problem presented in equation <ref> is an optimisation using the shear data γ, in fact it is the reduced shear (equation <ref>) that Glimpse uses to recover κ <cit.>. As an extension, the Glimpse algorithm can also perform the joint reconstruction with reduced shear and flexion, a third-order weak gravitational lensing effect <cit.> (although no flexion data are available for our galaxy shear catalogue). 

As the prior knowledge in this reconstruction relates to the quasi-spherical clustering of bound matter, enforced through a sparsity prior in starlet space, this method should better reconstruct the smaller scale non-Gaussian structure than the Wiener filter. 



§ DATA AND SIMULATIONS




 §.§ Dark Energy Survey Science Verification Data


The shear data are from the 139 deg^2 SPT-E field of the public DES SV data. This initial test data set was taken during an observing run before the official start of the full science survey. The galaxy catalogue comes from the SVA1 (Science Verification) Data Release[http://des.ncsa.illinois.edu]. Due to changes to the catalogues before final release (more galaxy shear measurements are now available to us), the catalogue used in this work is not identical to that used by <cit.>, even when the same data selections are made. All maps are therefore new, and slightly different to the previously published SV map.

The photometric redshifts from five optical filters (grizY) were estimated using the Bayesian Photometric Redshifts (BPZ) code (, , & ). The final median depth estimates are g∼24.0, r∼23.0, i∼23.0 and z∼22.4 (10-σ galaxy limiting magnitude). The “background galaxies”, the ones from which the shear is measured, are taken in the range 0.6 < z_mean < 1.2. The z_mean value for each galaxy is the mean of the posterior probability distribution function (PDF) estimated using the BPZ code. The PDF for each galaxy is very broad, giving a total stacked PDF of background galaxies that extends beyond the [0.6,1.2] redshift range, as can be seen in figure <ref>.

Using the ngmix shape catalogue, we apply a selection of sva1_flag= 0 & ngmix_flag = 0 to obtain galaxies with a well-measured shear. The  ngmix catalogue contains corrections to measurement bias, in the form of “sensitivities”, which can be applied to a weighted ensemble of hundreds or thousands of galaxies, but which cannot be applied per galaxy (which is not ideal for mass mapping). The structure of equation <ref> implies that a multiplicative shear bias would lead to a convergence amplitude bias. Under the assumption that multiplicative shear bias will not vary across the survey area, we correct all measured ellipticities by the same debiasing factor


    ϵ_obs, i =ϵ_measured, i×s̅^-1 ,


where i is a galaxy index and s̅ ( ≈ 0.82 ) is the mean sensitivity correction from all galaxies in our ngmix-selected catalogue. The total number of galaxies after the redshift and shape measurement selection is 1,628,663.

For the Kaiser-Squires reconstruction, the shear measurements are binned into angular pixels in a 256 × 256 map, with average pixel size of 4.11  arcmin, using a sinusoidal projection with a centre at RA=71.0 deg. This is similar to the 5 arcmin pixel scale of the original <cit.> map. The choice of central RA for Kaiser-Squires is to minimise the mask in the square projection, which is a large source of systematic error. For the Wiener filter, where the mask is taken into account, the shear measurements are also binned into angular pixels in a 256 × 256 map, but sinusoidally projected with a central RA=81.3 deg, to make the square maximally isotropic. The Glimpse algorithm does not bin the input shear measurements, but requires a pixel scale for the reconstruction, which we set as 3  arcmin using its gnomonic projection centred on RA=76.95 deg and DEC=-52.23 deg.

 

 







 §.§ redMaPPer Clusters
 

Groups and clusters of galaxies are expected to trace the highest density regions in the foreground. They are luminous objects that correspond to regions of highly non-linear growth, where the density field has deviated from Gaussianity.

The public redMaPPer cluster catalogue <cit.> used the redMaPPer algorithm to optically identify clusters and to estimate each cluster's richness, λ_RM. The richness is defined as the sum of the membership probabilities over all galaxies within a scale radius (chosen to minimise the scatter in the mass-richness relation); it gives an estimate for the number of galaxies in a cluster. Cluster mass is expected to scale approximately linearly with richness. The redshift uncertainty is excellent, around σ_z / (1+z) ∼ 0.01, due to the clusters containing large numbers of well modelled, red galaxies. The public redMaPPer catalogue used in this work contains only clusters with λ_RM⩾ 20, so that the clusters with less certainty of detection and characterisation are not used.

 

 



 §.§ Simulations


To compare the reconstructions between different methods, we use a simulated catalogue with a known true convergence. We use a set of N-body simulations developed for the DES collaboration and designed to be representative of the DES data <cit.>. The simulations used are N-body light cones composed from three boxes (1400^3, 2048^3, and 2048^3 particles in boxes of comoving length 1050 Mpc/h, 2600 Mpc/h, and 4000 Mpc/h respectively). The cosmological parameters for the simulations are: Ω_m = 0.286, Ω_Λ = 0.714, Ω_b = 0.047, σ_8 = 0.82, h_0 = 0.7, n_s = 0.96, w = -1. We apply a mask to match the SV data.

Source galaxies have randomly-assigned positions in the simulations, as correlation between the background galaxy positions and the weak lensing shear signal is expected to be negligible. The simulated catalogues contain the lensing matrix components, 𝒜_ij, for each galaxy, calculated with the ray-tracing code CALCLENS <cit.>. This provides the true κ and γ per galaxy, from which we derive the reduced shear. The shape noise due to the intrinsic ellipticities of the source galaxies, ϵ_s, is simulated by adding an ellipticity component to the reduced shear. Each noise realisation is generated from the data by randomly exchanging the ellipticity values between galaxies in the catalogue to remove the weak lensing signal and leave the shape noise.

We attempt to match the redshift distribution of the simulated galaxies to the observed redshift distribution, n(z). We use the stacked posterior probability density functions of individual galaxy redshifts from the selected data catalogue (figure <ref>), giving an estimate of the true underlying distribution. This assumes that


    n(z) = ∑_i p_i(z)  ,


where p_i(z) are the individual probability distributions for the galaxies from BPZ. This is not necessarily exact, due to errors in p_i (z) per galaxy <cit.>, but is a reasonable choice for a simulated catalogue. Using rejection sampling in bins of Δ z = 0.02 we select galaxies with a probability equal to the ratio between the desired n(z) from the data and the distribution in the simulation. One typical simulated catalogue contained 1,629,024 galaxies, slightly different to the data catalogue due to the sampling scheme, but with the desired n(z).



§ RESULTS


To ensure that the mass map tests are consistent with different output formats, all maps were converted onto a spherical pixelisation using HEALPix <cit.>. A HEALPix map comprises twelve subdivisions on the sphere, which are then each partitioned into NSIDE× NSIDE grids. Each pixel of a HEALPix supersampled NSIDE=4096 map was filled according to the value at the corresponding RA and DEC in the reconstructed maps. The supersampled high NSIDE maps were then degraded to NSIDE=1024. The true convergence maps from the simulations were directly binned from the convergence values at galaxy positions to NSIDE=1024. For all maps the same mask is applied, where pixels with no galaxies are masked.

Figure <ref> shows the mass map reconstructions from the SV shear data using the three different methods. An example simulation with truth and the three reconstructed maps is shown in figure <ref>. The “tuning parameters”, σ_smooth = 10.0  arcmin for Kaiser-Squires and λ = 3.0 for Glimpse, are tuned to maximise the Pearson correlation coefficient r with the underlying truth when tested on simulations.

Using a suite of 10 simulations, in Sec. <ref> we calculate the Pearson correlation coefficient between the truth and the reconstruction with different methods as a test of the reconstruction's quality. In Sec. <ref>, we calculate the root-mean-square error of the residuals between the truth and the reconstruction. In Sec. <ref> we calculate the variance of the 1-point distribution of the pixel values in the reconstruction and compare with the truth. In Sec. <ref> and Sec. <ref> we quantify the quality of the reconstruction of the phase and peak statistics respectively, by comparing to the simulated truth. The final result presented in Sec. <ref> compares the reconstruction from the DES SV shear data with foreground galaxy clusters from the redMaPPer catalogue (which are expected to trace non-linearities in the underlying density field). 

In this work we do not use correlation functions as a test of the map reconstruction. None of the mass mapping methods here are expected to reproduce the correct correlation functions or power spectra. It is simple to show this analytically with the Wiener filter, where despite the filter giving the MAP pixel values, the pixel variance, and therefore the power spectrum, is suppressed.

 



 §.§ Pixel Cross Correlation
 

We quantify the correlation between the true convergence from simulation and the reconstructed convergence of the simulated catalogue using the Pearson correlation coefficient. As with other metrics of success for mass map reconstruction, this can be used to tune the sparsity λ parameter and the smoothing scale for Kaiser-Squires.

The Pearson correlation coefficient, r, between the pixels' true convergence, κ^truth, and the reconstruction, κ^recon, is given by


    r =∑^n_i=1(κ_i^truth - κ̅^truth)(κ_i^recon - κ̅^recon)/√(∑^n_i=1(κ_i^truth - κ̅^truth)^2)√(∑^n_i=1(κ_i^recon - κ̅^recon)^2) ,


where the summations are over all pixels i in the map and κ̅ is the mean convergence in the map. 

In the left panels of figure <ref>, the Pearson r value from 10 simulations is plotted for varying tuning parameters. Almost all of the simulations and also their mean have a maximal Pearson r value at σ_smooth = 10.0  arcmin for Kaiser-Squires and at λ = 3.0 for Glimpse. 

Table <ref> presents the mean value from the 10 simulations, where the tuning parameter is chosen to maximise r when relevant. All methods show good correlation with the underlying true convergence. Both the Wiener filter and Glimpse have the same highest value of r = 0.37, 12 per cent higher than Kaiser-Squires. 

Note that the Pearson correlation coefficient as presented in equation <ref> is invariant under a rescaling of the reconstruction. Despite the Wiener filter reconstruction having values closer to zero, the Wiener filter maps still have good correlation to the truth. This second aspect is addressed in Sec. <ref> and in the second column of table <ref>.



 §.§ Pixel Residuals
 

The difference between the true convergence from simulation and the reconstruction in pixel i is defined as


    Δκ_i = κ_i^truth - κ_i^recon .


We define the root-mean-square error (RMSE) as


    RMSE(κ^truth, κ^recon) = √(1/n∑_i=1^n Δκ_i^2)


where n is the number of pixels. 

A smaller value of RMSE for a given method implies a better reconstruction according to this metric. It is this RMSE that the Wiener filter attempts to minimise using a linear filter, as defined in equation <ref>, by using an assumed signal covariance ⟨κκ^†⟩ (see Sec. <ref>).

The centre panel of figure <ref> shows that increasing the smoothing scale, σ_smooth, for Kaiser-Squires or the regularisation parameter, λ, for Glimpse initially reduces the pixel RMSE, but increased filtering contributes little beyond σ_smooth = 10.0  arcmin for Kaiser-Squires or λ = 3.0 for Glimpse. 

The smallest mean pixel RMSE is 1.0 × 10 ^-2 for Kaiser-Squires and 9.9 × 10 ^-3 for Glimpse. The Wiener filter, whose smoothing is constrained by the prior on C_ℓ and which therefore cannot be tuned, has a pixel RMSE of 9.4 × 10^-3.

 



 §.§ Pixel 1-Point Variance
 

The 1-point distribution can be thought as a histogram of the pixel values. Figure <ref> shows an example of such a histogram (derived from the simulated truth map and reconstructions of figure <ref>).

The mean of this distribution is unconstrained by weak lensing, due to an integration constant in equation <ref>. The variance of the 1-point distribution is increased compared to the underlying truth due to shape noise in the unsmoothed Kaiser-Squires reconstruction. A reconstruction method would aim to reduce the variance of the 1-point pixel distribution to match that of the underlying truth.

We define the estimate of the variance of the 1-point distributions of the truth or reconstructed κ as


    Var_truth =    1/n-1∑^n_i=1(κ_i^truth - κ̅^truth)^2 
    Var_recon =    1/n-1∑^n_i=1(κ_i^recon - κ̅^recon)^2  ,


where the notation matches equation <ref>. The ratio of these variances is given by


    Var_recon/Var_truth = ∑^n_i=1(κ_i^recon - κ̅^recon)^2 /∑^n_i=1(κ_i^truth - κ̅^truth)^2  ,


The closer this value is to 1, the better the variance of the pixel distribution matches the truth. Using 10 simulations we can calculate this quantity for different reconstruction methods (and at different smoothing scales or λ regularisation values where relevant). 

In figure <ref> the right panel shows the result of this test for Glimpse and Kaiser-Squires. Both methods show a pixel distribution that has too high variance for insufficient filtering, and too low variance for over-filtering. For Kaiser-Squires, the ratio is closest to 1 at a smoothing scale of σ_smooth = 5  arcmin. For Glimpse, the ratio is closest to 1 at a sparsity regularisation value of λ = 2.



Both of these reconstruction methods have a matching variance at a smoothing parameter value less than that which maximises the Pearson correlation coefficient r. If one chose this parameter to maximise the Pearson r value, such that λ = 3 and σ_smooth = 10  arcmin, a good reconstruction should also have the ratio of the variances as close to 1 as possible. 

The right column of table <ref> gives the mean variance ratio from 10 simulations with the different methods. The choice of λ = 3.0 and σ_smooth = 10  arcmin are the tuning parameters that maximise the Pearson r value for Glimpse and Kaiser-Squires respectively. Though Glimpse and the Wiener filter reconstructions both have the same Pearson r value, the variance of the pixel values of the Wiener filter is much lower with respect to the underlying truth than is the case for Glimpse. This can also be seen in the reconstructions of figure <ref>, where the Wiener filter pixel values are closer to zero than the simulated true convergence.

The histogram of figure <ref> shows, for one single example, the distributions matching what the results of the second column of table <ref> describe. Glimpse outperforms the other methods at matching the variance of the underlying truth, however it still falls short. Also, all methods, including Glimpse, have distributions which are symmetric, unlike the asymmetric, heavy-tailed distribution of the true κ values. 

Though Glimpse reconstructs maps with the 1-point distribution variance closest to the truth, it is also the only method to have convergence values dropping below the truth. These unphysical “negative peaks” can also be seen in the map reconstructions from data (figure <ref>) and from simulated catalogues (figure <ref>), and are likely to come from enforcement of sparsity for positive and negative wavelets equally. The physical motivation for Glimpse comes from a density field of superimposed halos. Though there should be no negative halos, negative wavelets are included to map the underdense regions, clearly at the expense of producing these very negative regions.



 §.§ Phase reconstruction
 

The summation over all m modes at each ℓ multipole in the angular power spectrum (equation <ref>) loses all phase information; only the magnitudes are retained. This phase information corresponds to the spatial distribution of anisotropies. As the phases are dependent on the physical underlying structure, they contain information beyond what can be gained by 2-point statistics. Their retention is a well-motivated, desired property of a mass mapping reconstruction. 

Inspired by <cit.>, who use phases to test the reconstruction after foreground removal from simulated Epoch of Reonization 21-cm maps, we use the phase residual as a metric of success between our three methods.

The phase difference between the true map and the reconstruction is defined as


    Δθ_ℓ m   = θ_ℓ m^truth - θ_ℓ m^recon
       = arg(a^truth_ℓ m) - arg(a^recon_ℓ m)  , 
         where  arg(z) = arctan(Im(z)/Re(z))  .


A small phase difference Δθ_ℓ m between the truth and the reconstruction implies that the phase has been well reconstructed. For random variables drawn from a Gaussian distribution, this would correspond to a small standard deviation. Here, however, a Gaussian distribution would be an inappropriate choice as it assumes the data are defined on an unbounded Euclidean space.

The two dimensional data space of phase pairs, {θ_ℓ m^truth, θ_ℓ m^recon}, is a torus, T^2, and the projected data space of the phase difference, Δθ_ℓ m, is a circle, S^1. On a circle, the maximum entropy, least informative, distribution for specified mean and variance is the von Mises <cit.>, which in one dimension is given by


    Pr(Δθ_ℓ m | C, μ) = 1/2  π  I_0 ( C )exp[ C cos (Δθ_ℓ m - μ)  ]  ,


where I_0 is the modified Bessel function of order 0, and C is a concentration parameter. For μ = 0, a large concentration parameter (analogous to 1/σ^2) would correspond to a small dispersion in the phase reconstruction error. The aim is therefore to compare the inferred value of the concentration, C, between different mass mapping methods, with a larger value of C implying a better phase reconstruction. 

By assuming that the error on the phase reconstruction is independent between phases, we can say that the phase differences, Δθ, are independent and identically distributed random variables, with a likelihood distribution given by


    Pr(Δθ | C, μ)    = ∏_ℓ m1/2  π  I_0 ( C )exp[  C cos (Δθ_ℓ m - μ)  ]
        = 1/[2  π  I_0 ( C )]^nexp[ C ∑_ℓ mcos (Δθ_ℓ m - μ)  ]  .


As only the relative values of C are needed to compare different mass mapping methods, the full posterior distribution is not required. Additionally, any reasonable prior distribution, Pr(C), will be either flat or monotonically decreasing above zero, so the ranking of maps by the largest maximum likelihood value or maximum posterior value of C will be identical. For the purposes of this comparison the simpler maximum likelihood estimate, Ĉ_MLE, will therefore do.

 

We calculate the maximum likelihood values of μ and C by taking the spherical harmonic transform of our HEALPix map to recover the a_ℓ m coefficients up to ℓ_max = 1024, calculating the phase residual as defined by equation <ref> between the truth and the reconstruction for each coefficient, and then maximising the likelihood (equation <ref>). The maximisation is performed using the scipy package BFGS algorithm (, , ), using 3 random initialisation values to test for robustness.

Figure <ref> show the results for the phase reconstruction from 10 simulations using Kaiser-Squires and Glimpse with varying tuning parameters. For Kaiser-Squires the mean phase reconstruction value, Ĉ_MLE, is maximised at σ_smooth=5.0 arcmin. For larger smoothing scales the phase reconstruction quality drops, as phase information is lost. For Glimpse the mean phase reconstruction value, Ĉ_MLE, is maximised at λ = 3.0. The maximum value of Ĉ_MLE is not particularly pronounced, and the Ĉ_MLE values are quite stable over a range of λ.



Table <ref> presents the mean values of Ĉ_MLE with the best tuning parameters for the three map reconstruction methods. Both Glimpse and the Wiener filter do much better than Kaiser-Squires for reconstructing the phases. Though the variance from these 10 different simulations is large, the Wiener filter does slightly better than Glimpse, as can be seen in figure <ref>.



 §.§ Peak Statistics
 

Peak statistics are a promising method for inferring cosmological parameters from data, as they access information beyond what can be inferred from 2-point correlation functions. Unlike higher order correlation functions, such as the bispectrum, peak statistics are inherently high signal-to-noise. They also probe the highly non-linear regions, where non-Gaussianity is greatest. The effect of masking is trivially taken into account by applying the identical mask to the suite of simulations used to construct a likelihood.

We cannot truly test which mass mapping method best constrains cosmology with the statistics of density peaks without fully deriving the posterior probability distributions of cosmological parameters. It is possible to test which method returns peaks which are distinguishable from noise and at which convergence values. Distinguishing a large number of peaks from noise at high values of κ would mean the map is reconstructing the non-linear regions well.

For a given convergence map, we can define a function, n(κ), that gives the number of peaks as a function of convergence. For a given mass reconstruction method we can compare the peaks in reconstructions from simulated data with the peaks in reconstructions from catalogues of “randoms”, with shape noise but no weak lensing shear signal (equivalent to γ = 0 in equation <ref>). If a given map from data or from a simulated catalogue has the same n(κ) as the random catalogues, then the mass mapping method used has been useless for peak statistics. On the other hand, if the map from data or simulation has a very different n(κ) function to that from the reconstruction from the random catalogue, then the map reconstruction method has recovered “true”, physical κ peaks.

 

 

In the DES SV cosmology constraints from peak statistics, <cit.> use this difference as the data vector used to constrain cosmology,


    n_Δ(κ_i) = n_data(κ_i) - n_randoms(κ_i)  .


This function is far from zero at a given κ if there is a large difference between the number of peaks counted in maps reconstructed (a) from data and (b) from random catalogues.

It is reasonable to believe that the number of peaks, n(κ_i), in the i^th bin, κ_i, is drawn from a Poisson distribution. The difference between two Poissonian random variables follows the Skellam distribution. Using this distribution, we expect the difference in the number peaks in maps from data and from random catalogues to have a mean given by


    μ_Δ(κ_i) = μ_data(κ_i) - μ_randoms(κ_i)  ,


and a variance given by


    σ^2_Δ(κ_i) = μ_data(κ_i) + μ_randoms(κ_i) .


We can therefore define a peak signal-to-noise estimate


    SNR(κ_i) = μ_Δ(κ_i)/√(σ^2_Δ(κ_i) )  .


Figure <ref> shows n_data, n_randoms, and n_Δ from Glimpse (λ = 3) from 10 simulations and 10 random catalogues. Here we define a peak as a local maxima in the HEALPIX map. Across different methods and smoothing parameters, the predicted variance from equation <ref> matches well with the estimated sample variance, verifying that the peak distribution is indeed Poissonian for a given κ.

Figure <ref> shows the peak signal-to-noise (SNR) estimates from 10 simulations and from 10 random catalogues as a function of κ and smoothing scale, for Kaiser-Squires, or λ, for Glimpse. As the peaks in the maps from data have higher convergence values than those from random catalogues, the SNR(κ) function is negative for low values of κ. 

In the figures, the Glimpse reconstruction gives better signal-to-noise estimates on the peaks than does the Kaiser-Squires reconstruction. For Kaiser-Squires, the largest positive and negative signal-to-noise values are 1.52 and -1.28. For Glimpse, the largest positive and negative are signal-to-noise values of 2.32 and -13.72. For the Wiener filter these values are 4.20 and -5.41.

The Wiener filter therefore has the highest signal-to-noise of the peak function n_Δ (κ), though the κ values of these peaks are very low. As can be seen in the reconstruction from the SV data in figure <ref>, the pixel values of the Wiener filter are much closer to zero. This is reflected in the peak statistic signal-to-noise values. In the left panel of figure <ref>, the Wiener filter detects negligibly few peaks with κ > 0.0125, whereas Glimpse detects peaks with positive signal-to-noise up to higher values of κ. It is at these high values where the non-Gaussian information due to non-linear structure formation can be probed.



 §.§ Foreground Clusters
 

Comparisons with foreground clusters of galaxies is an independent test of the mass map reconstructions, as it uses data (unlike our tests on simulations).

In figure <ref> the redMaPPer clusters described in Sec. <ref> are overlaid on the DES SV κ map reconstructions shown in figure <ref>. The maps show good spatial correlation between the locations of the clusters and the κ peaks in the map. 

The size of a cluster marker is the effective lensed cluster richness λ^eff_RM, rather than the redMaPPer cluster richness. This concept is adapted from the definition of κ_g presented in <cit.>. For a given cluster, this is defined as


    λ^eff_RM = p(z) ω(z)/a(z)×λ_RM×⟨λ_RM⟩/⟨λ^eff_RM⟩  ,



where z is the redshift of the cluster, p(z) is the lensing efficiency at the location of the cluster (see figure <ref>), and ω(z) is the comoving distance to the cluster (so that the first term matches the integrand of equation <ref>). The final term normalises the mean, where ⟨λ_RM⟩ is the average richness over all galaxy clusters. The effective lensed cluster richness gives the richness as “seen” by the lensing effect, where clusters at the peak of the lensing efficiency should contribute more to the κ map. We therefore calculate the correlation between λ^eff_RM for each cluster and the reconstructed κ value at the cluster centre. 

This method does not take into account multiple clusters overlapping in a given line of sight. In figure <ref>, many small clusters overlap on large peaks in the reconstructed κ map. The naive one-to-one correspondence between cluster and κ would mistake this for an excess of κ in the reconstruction. However, all methods will suffer equally from this assumption. A more thorough treatment of this overlapping effect is left for future work.



Table <ref> presents the Pearson correlation coefficient r between the λ^eff_RM value of each cluster and the κ^recon value at the corresponding pixel. The tuning parameters for Kaiser-Squires and Glimpse are chosen to maximise the Pearson correlation coefficient r between the reconstruction and the truth from simulations (see Sec. <ref>). 

Though both Glimpse and the Wiener filter take into account the noise and the mask in the data, and therefore do better than Kaiser-Squires, the Glimpse reconstructions show higher correlation with the effective richness of the foreground clusters than do the Wiener filter reconstructions. This is no surprise, as Glimpse is expected to do better at reconstructing non-Gaussian κ, which would correspond to the non-linear matter structures in which clusters of galaxies form.

 



§ CONCLUSIONS


In this work we have presented convergence map reconstructions using the public DES SV shear data with three different methods: Kaiser-Squires, Wiener filter, and Glimpse. Kaiser-Squires is a simple inversion from shear to convergence, whereas the Wiener filter and Glimpse use prior knowledge about the true convergence to help regularise the reconstruction and to reduce the effects of noise and missing data. The Wiener filter is a Bayesian MAP estimate if the signal and noise are Gaussian and the respective covariance matrices are known. The Glimpse method enforces a sparsity-promoting l_1 norm in a wavelet space where the wavelets represent positive, isotropic, quasi-spherical objects well. Glimpse is therefore expected to do well at reconstructing non-linear structures. The Wiener filter and Glimpse therefore aim to reconstruct different regimes: the linear and non-linear density field. 

The three methods were applied to realistic simulations of the DES SV shear data, for which an underlying true convergence is known. Using these simulations we are also able to tune the Kaiser-Squires smoothing scale, σ_smooth, and the Glimpse sparsity regularisation parameter, λ.

With these simulations we measure the Pearson correlation coefficient, r, between the truth and the reconstruction with different methods. Compared to the Kaiser-Squires reconstructions we find a 12 per cent improvement in Pearson correlation with both the Wiener filter and Glimpse. The tuning parameters of σ_smooth = 10  arcmin for Kaiser-Squires and λ = 3 for Glimpse maximise the Pearson correlation. We also measure the variance of the 1-point distribution of the reconstructed convergence. The Wiener filter suppresses the variance to 6.3  per cent of the truth, Kaiser-Squires to 37  per cent and Glimpse to 50  per cent of the truth. The tunable parameters here were those which maximised the Pearson correlation with the truth.

A large motivation for creating these maps is to reconstruct the convergence while still retaining the non-Gaussian information (which cannot be accessed with 2-point statistics such as the power spectrum). As such, we test the reconstruction of the harmonic phases, which is averaged out in the power spectrum, and the signal-to-noise of a peak statistic data vector, which is a popular probe of non-Gaussian information. The phase residuals between the truth and the reconstruction have the highest von Mises concentration with the Wiener filter (Ĉ_MLE = 0.591), with the Glimpse reconstruction performing comparably (Ĉ_MLE = 0.584). Both methods outperformed the Kaiser-Squires reconstruction (Ĉ_MLE = 0.501). 

With realistic data vectors for peak statistics generated from simulations, the maximum signal-to-noise value was increased by a factor of 3.5 for the Wiener filter and by a factor of 9 for Glimpse, compared to Kaiser-Squires. The signal-to-noise of the peak statistic data vector (n_Δ (κ)) is shown in figure <ref>, where Glimpse has significant signal-to-noise with high convergence peaks, where non-linearities in the underlying density field are highest. We predict these high value peaks are most useful for constraining cosmology beyond Gaussianity. In order to constrain cosmology with these different reconstruction methods, realistic simulations with different cosmological parameters or models must be used and the same reconstruction method should be applied to the simulations and data. As seen from our results, different reconstruction methods can produce convergence maps with different properties.

Finally, we switched from using simulations to instead using real observations (DES SV data). Here we measured the correlation between the reconstructed maps and the effective richness of the foreground redMaPPer clusters (this is the cluster richness as “seen” by the lensing effect). Table <ref> shows the results. Compared with Kaiser-Squires, the Wiener filter shows a 18 per cent increase and Glimpse shows a 32 per cent increase in correlation. This demonstrates with independent, cosmological data the ability of the methods to reconstruct non-linear structures.

The metrics we have used for comparing the three reconstruction methods are generic, and they have been inspired by recent applications of weak lensing mass maps to cosmological studies (e.g. , ). These metrics may not be optimal for evaluating every application of mass maps. Future studies can compare the efficiency of the three and other methods in end-to-end analyses; for example, with the estimation of cosmological parameters or identification of galaxy clusters.

Applying the Wiener filter and Glimpse methods to the DES Year 1 (Y1) shear catalogue would require extensions of the methods to account for the curved sky at large angular scales. The Y1 data covers ≈ 1500  deg^2 and contains ≈ 34,800,000 galaxies, so is a large increase in data volume from DES SV. This modification has already been done with an extension of Kaiser-Squires to the sphere by <cit.> for the Y1 DES data. These extensions would also be useful for the upcoming ≈ 5000  deg^2 DES Y3 shear catalogue.

Of future interest would be to use the Wiener filter or Glimpse convergence maps for scientific results, as we have shown that they reconstruct the convergence better than Kaiser-Squires according to many different metrics. 


 

We have made our map reconstructions (as shown in figure <ref>) available at <des.ncsa.illinois.edu/releases/sva1>.


 



§ ACKNOWLEDGEMENTS


NJ, FBA and J-LS acknowledge support from the European Community through the DEDALE grant (contract no. 665044) within the H2020 Framework Program of the European Commission. OL acknowledges support from a European Research Council Advanced Grant FP7/291329 and support from the UK Science and Technology Research Council (STFC) Grant No. ST/M001334/1. FBA also acknowledges the support of the Royal Society for a University Research Fellowship.

DJJ acknowledges the support of the National Science Foundation, award AST-1440254.

We are grateful for the extraordinary contributions of our CTIO colleagues and the DECam Construction, Commissioning and Science Verification teams in achieving the excellent instrument and telescope conditions that have made this work possible.  The success of this project also relies critically on the expertise and dedication of the DES Data Management group.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. 

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We  acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

We use the visualization software package SkyMapper[<https://github.com/pmelchior/skymapper>] for the map figures.







mnras
 










§ AUTHOR AFFILIATIONS


^1 Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

^2 Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown, 6140, South Africa

^3 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

^4 Laboratoire AIM, UMR CEA-CNRS-Paris 7, Irfu, SAp/SEDI, Service d'Astrophysique, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France

^5 Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

^6 Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

^7 Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

^8 Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, 85748 Garching, Germany

^9 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, 81679 München, Germany

^10 Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA

^11 Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

^12 Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA

^13 Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

^14 Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049 Madrid, Spain

^15 CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

^16 Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

^17 Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil

^18 Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil

^19 Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

^20 National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA

^21 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain

^22 Institut d'Estudis Espacials de Catalunya (IEEC), 08193 Barcelona, Spain

^23 Institute of Space Sciences (ICE, CSIC),  Campus UAB, Carrer de Can Magrans, s/n,  08193 Barcelona, Spain

^24 Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA

^25 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

^26 Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

^27 Department of Astronomy/Steward Observatory, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

^28 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

^29 Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

^30 Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

^31 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

^32 Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

^33 Department of Physics, The Ohio State University, Columbus, OH 43210, USA

^34 Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA

^35 Australian Astronomical Observatory, North Ryde, NSW 2113, Australia

^36 Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

^37 Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

^38 Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

^39 School of Physics and Astronomy, University of Southampton,  Southampton, SO17 1BJ, UK

^40 Department of Physics, Brandeis University, Waltham, MA 02453, USA

^41 Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil

^42 Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831





§ INDICATOR FUNCTION
 

We define the indicator function i_𝒞 (as used in equation <ref>) of a set 𝒞 as

    i_𝒞(x) = 
    	0     x ∈𝒞
    
    	+ ∞    .








