Weak lensing mass reconstruction using sparsity and a Gaussian random field
J.-L. Starck, K. E. Themelis, N. Jeffrey, A. Peel, F. Lanusse
AAstronomy & Astrophysics manuscript no. main © ESO 2021February 9, 2021
Weak lensing mass reconstruction using sparsity and a Gaussianrandom field
J.-L. Starck , K. E. Themelis , N. Je ff rey , , A. Peel , and F. Lanusse AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, F-91191 Gif-sur-Yvette, France Laboratoire de Physique de l’Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris,Paris, France Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny,1290 Versoix, SwitzerlandFebruary 9, 2021
ABSTRACT
Aims.
We introduce a novel approach to reconstruct dark matter mass maps from weak gravitational lensing measurements. Thecornerstone of the proposed method lies in a new modelling of the matter density field in the Universe as a mixture of two components:(1) a sparsity-based component that captures the non-Gaussian structure of the field, such as peaks or halos at di ff erent spatial scales;and (2) a Gaussian random field, which is known to well represent the linear characteristics of the field. Methods.
We propose an algorithm called MCALens which jointly estimates these two components. MCAlens is based on an alter-nating minimization incorporating both sparse recovery and a proximal iterative Wiener filtering.
Results.
Experimental results on simulated data show that the proposed method exhibits improved estimation accuracy compared tostate-of-the-art mass map reconstruction methods. (cid:135)
Key words.
Cosmology: observations – Gravitational lensing: weak – Methods: numerical – Techniques: image processing
1. Introduction
In recent years, there has been increasing interest in exploringthe two most dominant components of the universe, namely darkmatter and dark energy. To this end, large-scale imaging andspectroscopic surveys are currently under development, such asthe Euclid mission (Laureijs et al. 2011), the Rubin Observa-tory Legacy Survey of Space and Time (Abell et al. 2009), andthe Roman Space Telescope (Spergel et al. 2015), that will mapthe sky with unprecedented accuracy. A prominent cosmologicalprobe for these surveys is weak gravitational lensing.Weak gravitational lensing measures correlations in thesmall distortions of distant galaxies caused by the gravitationalpotential of massive structures along the line of sight. Its impacton distant source galaxies is twofold: the galaxy shapes are mag-nified by a convergence field, κ , while the galaxies’ ellipticitiesare perturbed from their underlying intrinsic value by a shearfield, γ . To contribute constraints on cosmological parametersand models, two-point correlation functions of the shear havebeen utilised with considerable success (Kilbinger et al. 2013;Alsing et al. 2016b). This type of correlation is su ffi cient to sta-tistically describe the Gaussian structures of the lensing field,such as those expected to be present either in the early universeor in large scales that are less a ff ected by gravitational collapse.To capture non-Gaussian structures such as those expected insmaller scales at later time, higher order moments of the shearneed to be employed (Munshi & Coles 2017).On the other hand, the convergence density field is not ob-served directly as a result of the mass-sheet degeneracy (Bartel-mann & Schneider 2001; Kilbinger 2015). Physically, the con-vergence field reveals the projected total matter density along the line of sight, weighted by a lensing kernel in the mid-distance be-tween the observer and the galaxy sources. The density field isinhomogeneous, encompassing Gaussian-type large-scale struc-tures, as well as non-Gaussian features, such as peaks. To shedlight on the way the convergence density field constrains cos-mology, peak statistics (e.g. Jain & Van Waerbeke 2000; Mar-ian et al. 2011; Lin & Kilbinger 2015; Liu & Haiman 2016;Peel et al. 2017b; Fluri et al. 2018; Li et al. 2019; Ajani et al.2020) and higher order correlation functions and moments, suchas the Minkowksi functionals (e.g. Kratochvil et al. 2012; Shi-rasaki et al. 2012; Petri et al. 2013), have been applied directlyon mass maps. It is therefore essential for mass mapping meth-ods to preserve both Gaussian and non-Gaussian features duringthe reconstruction process.Mass mapping methods solve an ill-posed problem due tothe irregular sampling of the lensing field and the low signal-to-noise ratio on small scales. A widely used algorithm to per-form mass mapping is the Kaiser-Squires method (Kaiser &Squires 1993), which is expressed as a simple linear opera-tor in Fourier space. However, it is inevitable for this estima-tor to su ff er from poor results, since it does not take specialcare of the noise or missing data. A di ff erent approach, moti-vated also by the Bayesian framework, is that of Wiener filter-ing (Wiener 1949). In this approach a Gaussian random field isassumed as a prior for the convergence map, which is respon-sible for inserting some bias that prevents our estimate fromover-fitting (Zaroubi et al. 1995). Moreover, a recently proposedstate-of-the-art method is the Gravitational Lensing Inversionand MaPping using Sparse Estimators (GLIMPSE2D) algorithm(Lanusse et al. 2016). GLIMPSE2D is a highly sophisticated al-gorithm that takes advantage of the sparse regularisation frame- Article number, page 1 of 14 a r X i v : . [ a s t r o - ph . C O ] F e b & A proofs: manuscript no. main work to solve the ill-posed linear inverse problem. GLIMPSE2Dis based on sparse representations (i.e. wavelets), and is thereforewell designed to recover piece-wise smooth features. An analyt-ical comparison between these three estimators is provided inJe ff rey et al. (2018a).In this paper we propose to bridge the gap between the sparseregularisation method of GLIMPSE2D and the Wiener filteringmethod by modelling the matter density field in the universe us-ing both linear and non-linear characteristics. Specifically, weassume that the density field is modelled as a mixture of twoterms: (1) a non-Gaussian term that adopts a sparse representa-tion in a selected wavelet basis (Starck et al. 2015), and (2) aGaussian term that is modelled using a Gaussian random field(Elsner, F. & Wandelt, B. D. 2013; Horowitz et al. 2019). Thenon-Gaussian signal component is able to capture the non-linearcharacteristics of the convergence field, such as peaks, while theGaussian component of the signal is responsible for capturingthe lower-frequency characteristics of the underlying field, suchas smooth variations. To our knowledge, this is the first time thatthis mixture modelling is proposed for mass map reconstruction.This paper is structured as follows. In Sect. 2 we intro-duce the formalism of weak gravitational lensing and describethe mass map reconstruction problem. To this end, we providea brief overview of the state-of-the-art algorithms of Kaiser-Squires, Wiener filtering, and GLIMPSE2D. We then present ourproposed mass mapping method in Sect. 3. The method is novelin the sense that it exploits both sparsity in the wavelet domainas well as a Gaussian random field model. Section 4 illustratesthe enhanced estimation performance of the proposed methodby providing experiments conducted on both simulated and realdata. Notation : We use ( · ) ∗ to denote the Hermitian transposeof a matrix or the adjoint operator of a transform. With (cid:107)·(cid:107) and (cid:107)·(cid:107) we denote the (cid:96) and (cid:96) norm respectively, ( (cid:107) x (cid:107) = (cid:80) Ni = | x i | , (cid:107) x (cid:107) = x T x ). The determinant of a matrix or the ab-solute value of a scalar is denoted by |·| , while diag( x ) stands fora diagonal matrix, that contains the elements of vector x on itsdiagonal. Finally, R N is the N -dimensional Euclidean space, denotes the zero vector, the all-ones vector, and I K is the K × K identity matrix.
2. Weak lensing mass mapping
Gravitational lensing describes the phenomenon where the lightemitted from distant galaxies is deflected as it passes through aforeground mass distribution. The lensing e ff ect causes the im-ages of distant galaxies to be distorted, with the distortion beingproportional to the size and shape of the projected matter distri-bution along the line of sight. Specifically, the mapping betweenthe source coordinates, β , and the lensed image coordinates, θ , isgiven by the lens equation (e.g. Kilbinger 2015), β = θ − ∇ ψ ( θ ) , (1)where ψ ( · ) defines the lensing potential that conceals the deflec-tion of light rays by the gravitational lens. Under the Born ap-proximation, which assumes that the lensing potential is weak enough, we may linearize the coordinate transformation ofEq. (1) by utilising the Jacobian A = ∂ β /∂ θ as, β i = A i j θ j , (2)where A i j = ∂β i /∂θ j are the coe ffi cients of the amplification ma-trix A , and we assume the Einstein summation convention. The symmetrical matrix A can also be parameterised in terms of the convergence , κ , and the shear , γ , as, A = (cid:34) − κ − γ − γ − γ − κ + γ (cid:35) . (3)The convergence can then be defined as a dimensionless quantitythat relates to the lensing potential through the Poisson equation, κ =
12 ( ∂ ∂ + ∂ ∂ ) ψ = ∇ ψ, (4)while the shear is mathematically expressed as a complex field,whose components also relate to ψ , γ =
12 ( ∂ ∂ − ∂ ∂ ) ψ and γ = ∂ ∂ ψ. (5)From Eq. (3), we see that the convergence causes an isotropicchange in the size of the source image, since it appears in the di-agonal of A . In comparison, the shear causes anisotropic changesto the image shapes. The convergence κ can also be interpretedvia Eq. (4) as a weighted projection of the mass density fieldbetween the observation and the source. Factoring out the term(1 − κ ) in Eq. (3) leaves the amplification matrix dependent onthe reduced shear, A = (1 − κ ) (cid:34) − g − g − g + g (cid:35) , (6)which is directly measured in lensing surveys and it is definedas g = γ/ (1 − κ ). In the weak lensing limit, where γ, κ (cid:28) g (cid:39) γ .In this paper we are interested in recovering the convergence κ from the reduced shear data. This is an ill-posed inverse prob-lem due to the finite sampling of the reduced shear over a lim-ited area of the survey and the presence of shape noise in themeasurements. In the following we review some of the state-of-the-art weak lensing mass reconstructing algorithms, namelythe Kaiser-Squires, the Wiener filtering, and the GLIMPSE2Dmethods. Kaiser-Squires:
A theoretical framework for reconstructingconvergence maps from the observable weak lensing shear in theFourier domain was proposed in Kaiser & Squires (1993). As wehave seen in Eqs. (4) and (5), the convergence and shear are bothexpressed as second order derivatives of the lensing potential.Their interrelation via the lensing potential ψ is expressed via atwo-dimensional convolution (Kaiser & Squires 1993), γ ( θ ) = π (cid:90) R d θ (cid:48) D ( θ − θ (cid:48) ) κ ( θ (cid:48) ) , (7)where D ( θ ) = − / ( θ − i θ ) . This convolution is equivalentlyexpressed in Fourier space as the element-wise multiplication,˜ γ ( k ) = π − ˜ D ( k )˜ κ ( k ) , (8)where k is the wavevector, and the Fourier transform of the ker-nel D ( θ ) is given by,˜ D ( k ) = π k − k + ik k k + k , (9)with k and k being the two frequency components of k . Dis-cretising Eq. 7 and adopting matrix notation, we may considerthat the observed shear γ is generated as a linear combination of Article number, page 2 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field the convolution matrix A and the unknown underlying field κ ,i.e., γ = A κ + n , (10)where n is the statistical uncertainty vector associated with thedata. Based on Eq. 8, A can be decomposed in Fourier space as A = F P F * , where F denotes the discrete Fourier transform, F * is its adjoint, and P is the diagonal operator that defines theconvergence field-shear relation in Fourier space, namely˜ γ = P ˜ κ = k − k k + ı k k k ˜ κ (11)where k = k + k and ˜ κ = F κ .Equation (11) corresponds to a discretised version of Eq. (8).As it stands, the Kaiser-Squires inversion of Eq. 11 su ff ers fromseveral drawbacks. First, it is not defined for k =
0, which stemsfrom the mass-sheet degeneracy (i.e. the mean value of the con-vergence field cannot be retrieved). Next, it is ill-posed, becausetypically the shear field is a discrete under-sampling of the un-derlying convergence field. Also, it does not take into accountmasked data. Nonetheless, the Kaiser-Squires estimate is stillused in practice due to its simplicity.
Wiener filtering:
The Wiener filter was introduced in the1940s (Wiener 1949), and it is the optimal linear filter in the min-imum mean square sense that provides a denoised version of thedesired signal. The Wiener filtered estimate of the convergencemap can be expressed using the linear equation κ G = W γ , (12)where the matrix W is given by, W = ( A Σ κ A * + Σ n ) − A * Σ κ , (13)and Σ κ and Σ n are respectively the pre-defined covariance matrixof the convergence and the noise. When the noise is not station-ary, the Wiener solution is not straightforward, and an iterativeapproach is required. This will be discussed in the next section.From a Bayesian perspective, Wiener filtering is equivalentto the maximum a posteriori (MAP) estimator given that the con-vergence is a zero-mean Gaussian signal with covariance Σ κ . In-deed, assuming that shear measurements are distorted by uncor-related Gaussian noise, the likelihood function shares the noiseproperties, that is, p ( γ | κ , Σ n ) = | π Σ n | − exp (cid:34) −
12 ( γ − A κ ) ∗ Σ n − ( γ − A κ ) (cid:35) , (14)while the distribution of the Gaussian random field can be writ-ten as p ( κ | Σ κ ) = | π Σ κ | − exp (cid:34) − κ ∗ Σ κ κ (cid:35) . (15)Given Eqs. (14) and (15) above, the MAP estimator is expressedas κ G = arg max κ { p ( κ | γ ) } ∝ arg max κ { p ( γ | κ , Σ n ) p ( κ | Σ κ ) }∝ arg max κ exp (cid:34) − κ ∗ Σ κ κ −
12 ( γ − A κ ) ∗ Σ n − ( γ − A κ ) (cid:35) , (16)leading to minimize κ G = arg min κ (cid:110) (cid:107) γ − A κ (cid:107) Σ n + (cid:107) κ (cid:107) Σ κ (cid:111) . (17) It is easy to see that the minimum is attained at κ G = W γ , where W is the Wiener filter given in Eq. (12).The implementation of the Wiener filter requires the inver-sion of the matrix W that includes both a signal covariance ma-trix component and the noise covariance matrix component.The unknown κ G is assumed to be a Gaussian randomfield with a covariance matrix diagonal in Fourier space, Σ κ = F * C κ F , where C κ is a diagonal matrix with diagonal valuesequal to the theoretical power spectrum P κ . When the noise isstationary, its covariance matrix is also diagonal with diagonalelements equal to the noise power spectrum P n . In this case, thefilter solution is obtained in Fourier space by ˜ κ G = ˜ W ˜ γ , wherethe Wiener filter is ˜ W = P κ P κ + P n . In practice the noise is generallynot stationary and depends on the number of shear measurementsin the area related to a given pixel of the shear field. Therefore thenoise covariance matrix is diagonal in pixel space, not in Fourierspace, and the Wiener solution becomes more problematic to de-rive, requiring either making a wrong assumption (i.e. that thenoise is stationary) or inverting a very large matrix. This rendersits inversion computationally complex and prone to numericalerrors. To circumvent this computationally intensive operation, aForward-Backward (FB) proximal iterative Wiener filtering wasproposed in Bobin et al. (2012) for Cosmic Microwave Back-ground spherical map denoising, exploiting the property that thesignal and noise covariance matrices are diagonal in pixel andFourier space, respectively. Eq. (17) comprises two separableterms, f ( κ ) = (cid:107) γ − A κ (cid:107) Σ n and f ( κ ) = (cid:107) κ (cid:107) Σ κ . (18)Following Forward-Backward methodology (Starck et al. 2015),we can solve Eq. (18) by designing an iterative fixed point algo-rithm as κ k + = prox µ f ( κ kG + ∇ f ( κ kG )) , (19)which is known to converge when µ < / (cid:107) A * Σ n − A (cid:107) , Com-bettes & Wajs (2005).Computing the proximal operator in Eq. (19), we end up withthe following iterative Wiener filtering algorithm, – Forward step: t = κ n + µ A * Σ − n ( γ − A κ n ) (20) – Backward step: κ n + = F ∗ (cid:18) P κ (cid:16) P η + P κ (cid:17) − (cid:19) F t , (21)where t is an auxiliary variable, µ = min( Σ n ), P η = µ I and κ =
0. This algorithm is free from matrix inversions, since both Σ n used in Eq. (20) and P κ used in Eq. (21) are diagonal ma-trices. A similar algorithm, exploiting transformations betweenpixel and harmonic space, was proposed by Elsner, F. & Wan-delt, B. D. (2013) using the messenger field framework. Suchmethods can also be adapted to e ffi ciently sample from the pos-terior distribution of the unknown mass map (Alsing et al. 2017;Je ff rey et al. 2018b). The Wiener approach recovers the Gaussiancomponent of the convergence field well, but it is far from op-timal at extracting the non-Gaussian information from the data,as peak-like structures are suppressed. This has motivated thedevelopment of sparse recovery methods based on wavelets. Sparse Recovery:
It has been shown that sparse recovery us-ing wavelets is a very e ffi cient way to reconstruct convergencemaps (Starck et al. 2006; Leonard et al. 2012; Lanusse et al. Article number, page 3 of 14 & A proofs: manuscript no. main (cid:96) -norm regularization in a wavelet-basedanalysis sparsity model. Sparsity in the Discrete Cosine Trans-form (DCT) domain was also proposed to fill the missing dataarea in the convergence map (Pires et al. 2009).The GLIMPSE algorithm avoids any binning or smoothingof the input data that could potentially cause loss of information.The primary cost function isarg min κ (cid:40) (cid:107) γ − T P F * κ (cid:107) Σ n + λ (cid:107) ω (cid:12) Φ ∗ κ (cid:107) + i R ( κ ) (cid:41) (22)where T is the nonuniform discrete Fourier transform (NDFT)matrix, P is defined in Eq. (11), λ is a sparsity regularizationparameter, ω is a weighting vector, Φ ∗ is the adjoint operator ofthe wavelet transform, and i R is an identity function that drivesthe imaginary part of the convergence to zero. This cost functionis also generalised in GLIMPSE2D in order to replace the shearby the reduced shear, or to incorporate the flexion information.It has been shown that GLIMPSE2D significantly outperformsWiener filtering for peaks recovery, while Wiener filtering doesbetter on the Gaussian map content (Je ff rey et al. 2018a). Deep Learning:
Deep learning techniques have recently beenproposed and seem very promising (Je ff rey et al. 2020). The in-put of the neural network is not the shear field directly, but ratherthe Wiener solution. We could imagine that the closer we areto the true solution with standard techniques such as Wiener orsparsity, the better deep learning can improve the solution. Openquestions related to deep learning remain to be answered, suchas its generalization to cosmologies not present in the trainingdata set, or the potential bias introduced by using a theoreticalpower spectrum in the Wiener solution serving as input of theneural network.
3. Modelling with sparsity and a Gaussian randomfield
We have seen in Sect. 2 two di ff erent models: the first mod-elling the convergence map as a Gaussian random field, leadingto good recovery of the large scales of the convergence map butsurpression peak structures, and the second assuming the con-vergence map is compressible in the wavelet domain (i.e. sparsemodelling). The sparse recovery is clearly complementary to theWiener solution, since it recovers peaks extremely well but theGaussian content poorly.To address these limitations, it seems natural to introducea novel modelling approach, where the convergence field κ isassumed to comprise two parts, a Gaussian and a non-Gaussian: κ = κ G + κ NG . (23)The non-Gaussian part of the signal κ NG is subject to a sparsedecomposition in a wavelet dictionary, while the component κ G is assumed to be inherently non-sparse and Gaussian.The Morphological Component Analysis (MCA) was pro-posed in Starck et al. (2004); Elad et al. (2005) to separate twocomponents mixed in a single image when these componentshave di ff erent morphological properties. This looks impossible,since we have two unknowns and one equation, but it was shownthat it is sometimes possible to extract these two componentsif we can exploit their morphological di ff erences. This requires having di ff erent penalisation functions C G and C NG on each ofthese two components, and we need to minimisemin κ G , κ NG (cid:110) (cid:107) γ − A ( κ G + κ NG ) (cid:107) Σ n + C G ( κ G ) + C NG ( κ NG ) (cid:111) . (24)MCA performs an alternating minimization scheme: – Estimate κ G assuming κ NG is known:min κ G (cid:110) (cid:107) ( γ − A κ NG ) − A κ G ) (cid:107) Σ n + C G ( κ G ) (cid:111) . (25) – Estimate κ NG assuming κ G is known:min κ NG (cid:110) (cid:107) ( γ − A κ G ) − A κ NG ) (cid:107) Σ n + C NG ( κ NG ) (cid:111) . (26)Examples of such decompositions can be seen on the MCAweb page . A range of MCA applications in astrophysics can befound in Starck et al. (2003); André et al. (2010); Möller et al.(2015); Bobin et al. (2016); Melchior et al. (2018); Joseph et al.(2019); Wagner-Carena et al. (2020). The Gaussian component κ G We use here the standard Wiener modeling where κ G is assumedto be a Gaussian random field: C G ( κ G ) = (cid:107) κ (cid:107) Σ κ , (27)and the solution of Eq. (25) is obtained using the iterative Wienerfiltering presented in the previous section. The non-Gaussian component
There are di ff erent ways to use a sparse model in the MCAframework. The most obvious would be to use standard (cid:96) or (cid:96) -norm regularisation in a wavelet-based sparsity model, as it isdone in the GLIMPSE2D algorithm. This would give: C NG ( κ NG ) = λ (cid:107) Φ ∗ κ NG (cid:107) p , (28)where p = Φ is the wavelet matrix, and λ is the regular-ization parameter (Lagrange multiplier).After implementing this approach, we found that largewavelet scales and Fourier low frequencies are relatively close,leading to di ffi culties in separating the information. We havetherefore investigated another approach, which involves first es-timating the set Ω of active coe ffi cients—i.e. the scales and po-sitions where wavelet coe ffi cients are above a given threshold—typically between 3 and 5 times the noise standard deviation rel-ative to each wavelet coe ffi cient. Ω can therefore be seen as amask in the wavelet domain, where Ω j , x = ffi -cient detected at scale j and position x , i.e. when | ( Φ ∗ A * γ ) j , x | >λσ j , x , and 0 otherwise. The noise σ j , x at scale j and position x can be determined using noise realizations as in the GLIMPSEalgorithm. An even faster approach is to detect the significantwavelet coe ffi cients on Φ ∗ A * Σ n − γ instead of Φ ∗ A * γ . Thenoise is therefore whitened, as the noise factor A * Σ n − is Gaus-sian with standard deviation equal to unity and with a uniformpower spectrum. We implemented both approaches giving simi-lar results, though the second is simpler.Once this wavelet mask Ω is estimated, we can estimate thenon Gaussian component κ NG bymin κ NG (cid:110) (cid:107) Ω (cid:12) Φ ∗ (( γ − A κ G ) − A κ NG )) (cid:107) + C NG ( κ NG ) (cid:111) . (29) http://jstarck.cosmostat.org/mca Article number, page 4 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field with C NG ( κ NG ) = i R ( κ NG ).This changes the original formalism since the data fidelityterm is now di ff erent, but it presents a very interesting advantage.Once Ω is fixed, the algorithm is almost linear and only the posi-tivity constraint remains. Therefore, we can easily derive a goodapproximation of the error map, just by propagating noise andrelaxing this positivity constraint. This will be further discussedin the following. Similarly to the GLIMPSE method, a positivityconstraint is applied on the non-Gaussian component κ NG . Peaksin κ can be on top of voids, and therefore have negative pixelvalues. As peaks are captured by the non-Gaussian component,they are positive by construction in κ NG , but the convergencemap κ = κ G + κ NG can still be negative at peaks positions. Largerare the non-Gaussianities, more we can expect MCAlens to im-prove over linear methods such as the Wiener one.The prior signal auto-correlation of the Gaussian componentis included within the signal covariance term of the Gaussiancomponent. We encode no explicit prior auto-correlation for thenon-Gaussian signal and no explicit prior cross-correlation be-tween the Gaussian and non-Gaussian component. Clearly, suchcorrelations exist but including their contribution in the prior inthis framework would be extremely di ffi cult theoretically and inpractice. However, such correlations will still appear in the finalreconstruction, driven by the correlation information in the data. Algorithm 1
MCAlens algorithm Input : Shear map γ , γ , signal and noise covariance Σ κ , Σ n ,and detection level λ . Initialize: κ NG(0) = κ G(0) = Ω = µ = min( Σ n ) , P η = µ I . Calculate wavelet coe ffi cients: α = Φ ∗ A * Σ n − γ . ∀ j , x , set Ω j , x = i f | ( α ) j , x | > λ . for n = , . . . , N max − do —————– Find κ NG ——————– Calculate the shear residual: γ r = γ − A ( κ G( n ) + κ NG( n ) ). Calculate the sparse residual: s r = A * Σ n − γ r . Calculate the sparse residual in the mask: s mr = Φ ( Ω (cid:12) ( Φ ∗ s r )). Get the new sparse component: S = κ NG( n ) + s mr . Positivity constraint: κ NG( n + = [ S ] + . —————- Find κ G ———————– Calculate the shear residual: γ r = γ − A ( κ G( n ) + κ NG( n + ). Forward step: t = κ G n + µ A * Σ − n γ r . Backward step: κ G n + = F ∗ (cid:18) P κ (cid:16) P η + P κ (cid:17) − (cid:19) F t . end for return (cid:16) κ G(N max ) , κ NG(N max ) (cid:17) . We solve the recovery problem of Eq. (23) using a two-stepoptimization procedure. First, a gradient descent step to min-imise Eq. (29) to recover the non-Gaussian component κ NG . Thisis followed by an iteration of the iterative Wiener filtering to min-imise Eq. (17). Details of the method are given in the Algo. 1.The number of scales N s used to compute the wavelettransform of an N x × N y image is automatically derived by N s = int (log(min( N x , N y ))). The λ parameter is a detection level,which was fixed for all our experiments with real and simulateddata to 5, i.e. 5 times the noise standard deviation, which is aconservative threshold which gave excellent results. Much attention has recently been given to the estimation of er-rors or uncertainties with mass map products. For linear meth-ods such as Wiener or Kaiser-Squires, it is easy to estimate thestandard deviation (or root-mean-square, RMS) per pixel, just bypropagating noise realisations using the same reconstruction fil-ters. The uncertainty per pixel does not, however, give the prob-ability that a clump in a reconstructed image is true or only dueto some noise fluctuations. Another approach, closer to certainscience cases with maps, estimates the significance of clumps.In Peel et al. (2017a), Monte Carlo simulations were used toaddress the significance of clumps. In Repetti et al. (2019), a hy-pothesis test called BUQO was proposed to do the same task,requiring the user to define manually a mask around the clump.Similarly, in Price et al. (2019), hypothesis tests of Abell-520cluster structures using Highest Posterior Density Regions wereperformed.With MCAlens, we can include aspects of both approaches.
RMS and SNR maps
In Algo. 1, there are two steps involving a non-linear operator,first for the estimation of Ω in line 4 and the second in line 11to perform the positivity constraint. To propagate noise realiza-tions, Ω has to be set to the one obtained with the data, and thisnon-linearity step therefore does not occur, so only the positiv-ity remains. A full linear algorithm could therefore be obtainedjust by removing this positivity constraint during the noise prop-agation, by replacing κ NG( n + = [ S ] + in Algo. 1 line 11 by κ NG( n + = S . This leads to more noise entering the solution,since few pixel values with negative values in κ NG are not thresh-olded, and the derived RMS map is therefore slightly conserva-tive. Hence, we can build noise realizations, run the MCAlensalgorithm on each realization, and derive the RMS map by tak-ing pixel per pixel the standard deviation of the obtained recon-structed maps. The SNR map is derived by dividing the absolutevalue of the reconstructed map by the RMS map. Significance map
In Algo. 1 line 4, the wavelet mask Ω is obtained by comparingthe wavelet coe ffi cients, α = Φ ∗ A * Σ n − γ , to the threshold λσ where σ = H that the waveletcoe ffi cient is due to noise only, and if a given wavelet coe ffi cient α j , x is such that | α j , x | > λ , then the H hypothesis is rejected andwe assume that the coe ffi cient amplitude cannot be explained bynoise fluctuations and is therefore due to signal. λ is thereforedirectly related to the significance of the wavelet coe ffi cients,and the mask Ω indicates which coe ffi cients are detected witha given significance level. Ω j , x is binary, and we build the sig-nificance map s by s x = (cid:80) j Ω j , x , i.e. a simple coadding of allbinary scales. Such a map could also be used as a way to au-tomatically derive the user mask required in the BUQO method(Repetti et al. 2019). We present in Sect. 4 examples of RMS,SNR and significance maps. When a wide-field map needs to be reconstructed, the flat ap-proximation cannot be used anymore, and we have to build amap on the sphere. A traditional approach is to decompose the
Article number, page 5 of 14 & A proofs: manuscript no. main sphere into overlapping patches, assume a flat approximation oneach individual patch, reconstruct each patch independently, andfinally recombine all patches on the sphere. This solution is cer-tainly good enough to recover clumps relative to clusters, i.e.the non-Gaussian component, but certainly not for the Gaussiancomponent which contains information at low frequencies. Inthe framework of the DES project, a 1500 deg map has beenreconstructed Chang et al. (2018), using HEALPIX pixelisa-tion (Górski et al. 2005) and a straightforward spherical Kaiser-Squires algorithm consisting in:1. calculating a spin transform of the HEALPIX shear map toget both the E and B spherical harmonic coe ffi cients,2. smoothing by a Gaussian both E and B modes in the spheri-cal harmonic domain,3. applying an inverse spherical transform independently oneach of these two modes to get the two convergence maps κ E and κ B .Sparsity in a Bayesian framework (Price et al. 2020b) and for-ward fitting in harmonic space (Mawdsley et al. 2020) have alsobeen recently proposed to do spherical mass mapping. Usingsimilarly a HEALPIX shear and convergence pixelisation, wecan also easily derive an extension of both the iterative-Wienerfiltering and the MCAlens algorithm by – replacing the matrix A = F P F * in Eq. (10) by A = Y Y ∗ ,where s Y and s Y ∗ represent the forward and inverse spin-sspherical harmonic transforms respectively, – replacing the the wavelet decomposition Φ used in Eq. (28)by the spherical wavelet decomposition (Starck et al. 2006).Then the same MCAlens algorithm given in Algo. 1 can be usedto derive a spherical convergence map. An example is given inSect. 4.3. We focused earlier on the convergence map (i.e. E-mode). TheMCAlens algorithm given in Algo. 1 remains valid if one wantsto estimate both E and B mode by adopting the following nota-tion: γ = (cid:32) γ γ (cid:33) , κ = (cid:32) κ E κ B (cid:33) , A = k − k k k k k k k k − k − k k , α = (cid:32) α E α B (cid:33) , P κ = (cid:32) P κ E P κ B (cid:33) , Φ = (cid:32) Φ E Φ B (cid:33) ,where the matrix Φ consists in applying a sparse decomposi-tion independently on each mode. In practice we use the samewavelet decomposition for both (i.e. Φ E = Φ B ). The delicatepoint is the Wiener filter to apply to the B-mode at line 15 of thealgorithm. In theory, P κ B =
0, and by construction, no Gaussiancomponent can be recovered. Since the B-mode is mainly use-ful for investigation of systematic errors, we are not interestedin recovering the B-mode least square estimator (which is zero),and we find it more useful to process the B-mode similarly tothe E-mode. We therefore advocate rather to use P κ B = P κ E . Thisway the fluctuations in the E-mode can properly be compared tothose in the B-mode.
4. Experimental results
We use a simulation derived from RAMSES N-body cosmo-logical simulations (Teyssier 2002), with a Λ CDM model (see
Fig. 1.
RAMSES simulations: Error versus scale for Wiener in red andMCAlens in blue. ), with a pixel size of 0 . (cid:48) x 0 . (cid:48) and a galaxy redshift of 1. To get a realistic mask andnoise behavior, we use the MICE pixel noise covariance derivedfor the DES project (see Appendix A). As the pixel resolutionis di ff erent, it leads to an optimistic realization, but it has theadvantage of illustrating well the impact of our MCA model.We have run MCAlens using 100 iterations, λ = σ ), the MICE covariance matrix, and the used the-oretical power spectrum was the true map power spectrum. Toevaluate the results, we ran 100 di ff erent noise realizations, andwe applied both Wiener filtering and MCAlens. We calculatedthe reconstruction error at di ff erent resolutions with Err % ( σ ) = (cid:107) G σ ( M ( κ − κ t )) (cid:107)(cid:107) M κ t (cid:107) (30)where (cid:107) x (cid:107) = (cid:113)(cid:80) i x i , κ is the reconstruct convergence map, κ t the true convergence map, G σ ( x ) is the convolution of x with aGaussian with standard deviation σ , and M is the binary maskwith M k = k (i.e. we have data at this location) and 0 otherwise. Figure 1shows the mean error Err % ( σ ) for both the Wiener and MCAlenssolutions. The black curve shows the di ff erence between theWiener error and the MCAlens error, allowing to better visualisethat the improvement is larger at fine scales. It is interesting tonote that at MCAlens is better than Wiener even at large scales.Indeed, when the MCAlens non-Gaussian component contains asignificant amount of features as in the case of this experiment,these features have also a non negligible contribution on largerscales, which explain why MCAlens remains better than Wienerat large scales. MCAlens leads to a clear improvement in termsof quadratic error. The top panels of Fig. 2 show respectivelythe simulated convergence map and the MCAlens reconstructedmap. The bottoms panels show the Gaussian and the non Gaus-sian part recovered from the noisy the data. The sum of thesetwo components is equal to the MCAlens reconstruction (topright). This shows that both the non-linear and the linear com-ponents can be recovered well. Figure 3 shows respectively theRMS map, the SNR map and the significance map. Article number, page 6 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field
Fig. 2.
RAMSES simulations: Top, true convergence map and MCAlens recovery; bottom, Gaussian components and sparse components. The sumof these two maps is equal to the top right MCAlens map.
We use here a convergence map, released by theColumbia Lensing group (Liu et al. 2018), corre-sponding to a cosmological model with parameters (cid:110) M ν , Ω m , A s , M ν , h , w (cid:111) = { . , . , . , , . , − } , and witha pixel size of 0 . (cid:48) x 0 . (cid:48) . We rebinned the map to 0 . (cid:48) x0 . (cid:48) , and similarly to the previous experiment, we simulatednoisy data using the same MICE covariance, applying a globalrescaling in order to have realistic noise corresponding to amean number of galaxies equal to 30 per arcmin , as we expectin future space lensing surveys. To evaluate the results, we ran100 di ff erent noise realizations, and we applied Kaiser-Squires,sparse recovery, Wiener filtering, and MCAlens. For the sparserecovery and MCAlens, we used λ = σ ),and for Wiener and MCAlens we the used the theoretical powerspectrum as the true convergence map power spectrum.Figure 4 shows the error computed using Eq. (30). We cansee that, in addition to the well-recovered peaks, MCAlensalso leads to a clear improvement in terms of quadratic er- http://columbialensing.org ror compared the other methods. On the contrary to the RAM-SES experiment, we see here a convergence at large scales be-tween MCAlens and Wiener. This is due that the MCAlens non-Gaussian component contains only a few peaks (see Fig. 5, mid-dle right), which have therefore a negligible contribution to thelargest scales.The top row of Fig. 5 shows respectively the simulated con-vergence map, the Wiener result, and the Kaiser-Squires re-construction with smoothing applied. The middle shows theMCAlens map and its Gaussian and non Gaussian parts recov-ered from the noisy the data. The sum of these two componentsis equal to the MCAlens. Similarly to previous experiments, boththe non-linear and the linear components are well recovered. Thebottom row of Fig. 5 shows respectively the RMS map, the SNRmap and the significance map. We have created a full shear map from the full sky MICE simu-lated map and its covariance matrix, which has around 3 galaxiesper arcmin . We ran the spherical MCAlens method with 200 it- Article number, page 7 of 14 & A proofs: manuscript no. main
Fig. 3.
RAMSES simulations: from left to right, RMS map, SNR map and significance map.
Fig. 4.
Columbia simulations: error versus scale for four di ff erent meth-ods: Kaiser-Squires (light blue), Wiener (red), sparsity (orange) andMCAlens (blue). erations, λ = σ ), and we used the powerspectrum of the simulation as the theoretical power spectrum.The top row of Fig. 6 shows the simulated noise-free map at 1and 3 degrees, while the bottom row shows the MCAlens non-Gaussian and Gaussian components. In this last section, we apply MCAlens to reconstruct a conver-gence map of the 1.64 deg HST / ACS COSMOS survey (Scov-ille et al. 2007). In this work, we make use of the bright galaxiesshape catalogue produced for (Schrabback et al. 2010).The results after applying MCAlens on COSMOS data arepresented in Fig. 7. In the top row are the galaxy count map,the Wiener map, and the Kaiser-Squires map smoothed with aGaussian at a Full Width at Half Maximum of 2.4 arcmin. Thebottom row shows respectively the Glimpse, MCAlens E-modeand MCAlens B-mode maps. White dots show the locationsand redshifts of X-ray selected massive galaxy clusters from theXMM-Newton Wide Field Survey (Finoguenov et al. 2007) with0 . < z < .
5. Conclusion
A novel mass mapping algorithm has been presented that is ableto recover high resolution convergence maps from weak grav-itational lensing measurements. Our proposed process involvesa model with two components, a Gaussian and a non-Gaussian,for which we have developed an e ffi cient algorithm to derive thesolution. We have shown that we can also handle a non-diagonalcovariance matrix. We have extended the method so it can dealwith spherical maps, which is needed for future surveys such asthe Euclid space mission. Our experiments clearly show a sig-nificant improvement compared to the state of art.In the spirit of reproducible research, the MCAlens algo-rithm is publicly available in the CosmoStat’s Github package ,including the material needed to reproduce the simulated experi-ences (folder examples / mcalens_paper and script make_fig.py). Acknowledgment
We thank Tim Schrabbaack for sending us his COSMOSshear catalog and the Columbia Lensing group ( http://columbialensing.org ) for making their suite of simulatedmaps available, and NSF for supporting the creation of thosemaps through grant AST-1210877 and XSEDE allocation AST-140041. https://github.com/CosmoStat/cosmostat Article number, page 8 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field
Fig. 5.
Columbia convergence map recovery: Top, true convergence map, Wiener map, and Kaiser-Squires map smoothed with a Gaussian havinga Full Width at Half Maximum of 3.8 arcmin. Middle, MCAlens map and its Gaussian and sparse and. The sum of these two last maps is equal tothe first one. Bottom, RMS, SNR and significance maps. Article number, page 9 of 14 & A proofs: manuscript no. main
Fig. 6.
MICE simulations: top, true convergence map at 1 and 3 degrees. Bottom, MCAlens sparse and Gaussian components.
Fig. 7.
COSMOS data: Top, galaxies count map, Wiener map, and Kaiser-Squires map smoothed with a Gaussian having a Full Width at HalfMaximum of 2.4 arcmin. Bottom, Glimpse, MCAlens and MCAlens B-mode map.Article number, page 10 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field
Appendix A: MICE simulations
We use the public MICE (v2) simulated galaxy catalogue, whichis constructed from a lightcone N-body dark matter simula-tion (Fosalba et al. 2015; Crocce et al. 2015; Fosalba et al.2015; Carretero et al. 2015; Ho ff mann et al. 2015; Tallada et al.2020). The MICE catalogue provides the calculated weak lens-ing (noise-free) observables: shear and convergence. In a givenpatch of simulated sky we select galaxies in the redshift z range[0 . , . ∼ . Redshift due to the expansion of the Universe is used as a proxy fordistance to a galaxy, as it is an observable that monotonically increaseswith distance from an observer.
Uncorrelated, complex shape noise values are randomlydrawn from a Gaussian distribution and added to the shear valueof each selected galaxy. This noise per galaxy is zero mean andhas variance 2 σ (cid:15) = . ff rey et al.2018a)). The final pixelised noise ( n ) has variance that dependson the number of galaxies per pixel.In our simulated data, we mimic these conditions by choos-ing to remove all galaxies in given regions. Here there are noshear measurements available and the noise variance is e ff ec-tively infinite.The top row of Fig. A.1 shows a simulated convergence mapwith the DES SV footprint and mask, and the bottom row showstwo simulated observed shear component maps, γ and γ . Appendix B: A Wiener tour
Appendix B.1: Wiener and inpainting
Missing data is a common problem for galaxy surveys as fore-ground objects obscuring the background galaxies have to be“masked out”. In addition to the non stationary noise, observedshear fields therefore also present missing data. Noting the mask M as equal to 1 if we have shear measurements at a pixel posi-tion and zero otherwise, missing data can be handled in Wienerfiltering by forcing the noise covariance matrix to be very highat locations where M =
0. This leads the solution to be di ff er-ent from zero and smoothed in the missing data area, and re-move border e ff ect artefacts. The Wiener filtering can thereforebe seen as an inpainting technique, since it fills the missing areain the image. Alternative inpainting techniques were proposed inthe past, through sparse recovery techniques (Pires et al. 2009;Lanusse et al. 2016) or Gaussian constraint realizations (Zaroubiet al. 1995; Je ff rey et al. 2018a). Since the Wiener model as-sumed the solution to be a Gaussian random field, it would makesense to have a solution where the inpainted area presents thesame statistical properties as in non inpainted aera. This prop-erty is by construction verified with constraint realizations, andit was shown it is also the case with sparse inpainting (Pires et al.2009). A similar sparse inpainting can be very easily included inthe proximal Wiener filtering, minimizing the following equa-tion: κ G = arg min κ (cid:110) (cid:107) M ( γ − A κ ) (cid:107) Σ n + β (cid:107) κ (cid:107) Σ κ + λ (cid:107) Φ ∗ κ (cid:107) p (cid:111) , (B.1)where p = Φ is the Discrete Cosine Dictionary (DCT),and λ is a Lagrangian parameter. We end up with the forward-backward algorithm, called InpWiener : – Forward step: t = κ n + µ A * Σ − n ( M ( γ − A κ n )) (B.2) – Backward step: κ n + = M u + (1 − M ) ∆ Φ ,λ u (B.3)where and ∆ Φ ,λ is the proximal operator defined in Pireset al. (2009) which consists in applying a DCT transform to u , threshold the DCT coe ffi cients and reconstruct an imagefrom the thresholded coe ffi cients, and u = F ∗ (cid:18) P κ (cid:16) P η + P κ (cid:17) − (cid:19) F t , (B.4)Areas where we have information are processed as in the usualWiener case, while the inpainting regularization impacts areawith missing data (i.e. when M = Article number, page 11 of 14 & A proofs: manuscript no. main
Fig. A.1.
Top, true MICE convergence map with DES SV footprint and mask. Bottom, Two simulated observed shear components, γ and γ . – KS: if β = , λ = Σ n is diagonal with constant val-ues along the diagonal (i.e. stationary Gaussian noise), In-pWiener leads to the non-iterative standard Kaiser-Squiressolution. – GKS: If β = λ =
0, the least square estima-tor is derived with the iterative algorithm: κ n + = κ n + µ A * Σ − n ( M ( γ − A κ n )), with µ = min( Σ n ). As it general-izes the Kaiser-Squires method, we will call this algorithm GKS . – FASTLens: If β = Σ n is diagonal with constant val-ues along the diagonal, InpWiener leads to the FASTLensinpainting algorithm (Pires et al. 2009). – GIKS: If β = InpWiener leads to an inpainted general-ized the Kaiser-Squires solution where the
InpWiener
For-ward is unchanged, and the Backward step becomes: κ n + = M t + (1 − M ) ∆ Φ ,λ t (B.5)Similarly to sections 3.5 and 3.5, these algorithms can handledata on the sphere and reconstruct jointly E and B modes. Inpainted Wiener Experiment
To test
InpWiener , we use the public MICE (v2) simulatedgalaxy catalogue presented Appendix A.Figure B.1 shows the Wiener solution (left) and the in-painted Wiener solution (right) derived from the shear compo-nents shown in Fig. A.1.
Appendix B.2: Agnostic Wiener Filtering
The Wiener method needs to know the theoretical power spec-trum P κ , and the solution therefore varies with the assumed cos-mological model used to derive P κ . To avoid this issue, a solu-tion could to estimate P κ directly from the shear measurements,for instance using a mask correction as in Upham et al. (2020).Bayesian techniques have also been used to infer both the mapand power spectrum (Wandelt et al. 2004; Jasche & Lavaux2014; Alsing et al. 2016a). An alternative approach is to usedthe GIKS inpainting algorithm to fill first the missing area andthen to compute the power spectrum of the inpainted map. Ap-plying
GIKS on the data and a set of R noise realizations, thefinal estimator is P κ = powspec( κ Data ) − R (cid:88) i powspec( κ Rea i ) . (B.6)Since the data noise P κ will still be noisy, a final denoising stepor a function fitting can be done.As an illustration, we fitted the function f ( k , a , u , e , c ) = exp( | a ∗ ( uk ) − e | ) + c to the estimated noisy P κ . Figure B.2 showsan example of an estimated power spectrum following this ap-proach. Article number, page 12 of 14.-L. Starck et al.: Weak lensing mass reconstruction using sparsity and a Gaussian random field
Fig. B.1.
Wiener and inpainted Wiener solutions.
Fig. B.2.
Theoretical power spectrum of one convergence map. In blackthe true theoretical power spectrum; in red, the spectrum estimated us-ing the
GIKS algorithm, corrected from noise power spectrum, and inblue the fit.
Experiment: Impact of an unknown theoretical power spec-trum
In this experiment, we used the same public MICE simulations,and we extracted 18 shear di ff erent shear maps with di ff erentnoise realizations. For each of them, we applied the forward-backward Wiener algorithm using the true theoretical powerspectrum, and we applied the inpainted agnostic Wiener methodon the same data, re-estimating for each of the 18 shear data setsthe theoretical power spectrum. Fig. B.3 shows the reconstruc-tion errors at di ff erent resolutions, for Kaiser-Squires, Wienerand inpainted agnostic Wiener. We can easily see that the in-painting has no impact on the reconstruction error, which is ex-pected since only the area where the mask is equal to one is usedin the error calculation, and also that the agnostic approach alsohas very little impact on the final solution. We do not claim thatwe should use an agnostic approach when applying a Wiener fil-tering, but it is interesting to have this option available, i.e. whenusing Wiener without any assumption about the cosmology, andto give us the possibility to check if cosmological priors impactthe results. Fig. B.3.
Reconstruction error at di ff erent resolution for Kaiser-Squires,Wiener and inpainted agnostic Wiener. References
Abell, P. A., Burke, D. L., Hamuy, M., et al. 2009, Lsst science book, version2.0, Tech. rep.Ajani, V., Peel, A., Pettorino, V., et al. 2020, arXiv e-prints, arXiv:2001.10993Alsing, J., Heavens, A., & Ja ff e, A. H. 2016a, Monthly Notices of the RoyalAstronomical Society, 466, 3272Alsing, J., Heavens, A., & Ja ff e, A. H. 2017, MNRAS, 466, 3272Alsing, J., Heavens, A., Ja ff e, A. H., et al. 2016b, MNRAS, 455, 4452André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102Bartelmann, M. & Schneider, P. 2001, Physics Reports, 340, 291Bobin, J., Starck, J. L., Sureau, F., & Fadili, J. 2012, Advances in Astronomy,2012, 703217Bobin, J., Sureau, F., & Starck, J. L. 2016, A&A, 591, A50Carretero, J., Castander, F. J., Gaztañaga, E., Crocce, M., & Fosalba, P. 2015,MNRAS, 447, 646Chang, C., Pujol, A., Mawdsley, B., et al. 2018, MNRAS, 475, 3165Combettes, P. L. & Wajs, V. R. 2005, Multiscale Modeling & Simulation, 4, 1168Crocce, M., Castander, F. J., Gaztanaga, E., Fosalba, P., & Carretero, J. 2015,MNRAS, 453, 1513Elad, M., Starck, J.-L., Donoho, D., & Querre, P. 2005, Applied and Computa-tional Harmonic Analysis, 19, 340–358Elsner, F. & Wandelt, B. D. 2013, A&A, 549, A111Finoguenov, A., Guzzo, L., Hasinger, G., et al. 2007, ApJS, 172, 182Fluri, J., Kacprzak, T., Sgier, R., Refregier, A., & Amara, A. 2018, J. CosmologyAstropart. Phys., 2018, 051Fosalba, P., Crocce, M., Gaztanaga, E., & Castander, F. J. 2015, MNRAS, 448,2987 Article number, page 13 of 14 & A proofs: manuscript no. main
Fosalba, P., Gaztanaga, E., Castander, F. J., & Crocce, M. 2015, MNRAS, 447,1319Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, Astrophysical Journal, 622[ arXiv:astro-ph/0409513 ], 759–771Ho ff mann, K., Bel, J., Gaztanaga, E., et al. 2015, MNRAS, 447, 1724Horowitz, B., Seljak, U., & Aslanyan, G. 2019, J. Cosmology Astropart. Phys.,2019, 035Jain, B. & Van Waerbeke, L. 2000, ApJ, 530, L1Jasche, J. & Lavaux, G. 2014, Monthly Notices of the Royal Astronomical Soci-ety, 447, 1204Je ff rey, N., Abdalla, F. B., Lahav, O., et al. 2018a, MNRAS, 479, 2871Je ff rey, N., Heavens, A. F., & Fortio, P. D. 2018b, Astronomy and Computing,25, 230Je ff rey, N., Lanusse, F., Lahav, O., & Starck, J.-L. 2020, MNRAS, 492, 5023Joseph, R., Courbin, F., Starck, J. L., & Birrer, S. 2019, A&A, 623, A14Kaiser, N. & Squires, G. 1993, ApJ, 404, 441-450Kilbinger, M. 2015, Reports on Progress in Physics, 78, 086901Kilbinger, M., Fu, L., Heymans, C., et al. 2013, MNRAS, 430, 2200Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2012, Phys. Rev. D, 85, 103513Lanusse, F., Starck, J.-L., Leonard, A., & Pires, S. 2016, A&A, 591, A2Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv preprint arXiv:1110.3193Leonard, A., Dupé, F.-X., & Starck, J.-L. 2012, A&A, 539, A85Li, Z., Liu, J., Matilla, J. M. Z., & Coulton, W. R. 2019, Phys. Rev. D, 99, 063527Lin, C.-A. & Kilbinger, M. 2015, A&A, 576, A24Liu, J., Bird, S., Matilla, J. M. Z., et al. 2018, J. Cosmology Astropart. Phys.,2018, 049Liu, J. & Haiman, Z. 2016, Phys. Rev. D, 94, 043533Marian, L., Hilbert, S., Smith, R. E., Schneider, P., & Desjacques, V. 2011, ApJ,728, L13Mawdsley, B., Bacon, D., Chang, C., et al. 2020, MNRAS, 493, 5662Melchior, P., Moolekamp, F., Jerdee, M., et al. 2018, Astronomy and Computing,24, 129Möller, A., Ruhlmann-Kleider, V., Lanusse, F., et al. 2015, J. Cosmology As-tropart. Phys., 2015, 041Munshi, D. & Coles, P. 2017, Journal of Cosmology and Astroparticle Physics,2017, 010Peel, A., Lanusse, F., & Starck, J.-L. 2017a, The Astrophysical Journal, 847, 23Peel, A., Lin, C.-A., Lanusse, F., et al. 2017b, A&A, 599, A79Petri, A., Haiman, Z., Hui, L., May, M., & Kratochvil, J. M. 2013, Phys. Rev. D,88, 123002Pires, S., Starck, J. L., Amara, A., et al. 2009, MNRAS, 395, 1265Price, M. A., Cai, X., McEwen, J. D., et al. 2020a, MNRAS, 492, 394Price, M. A., McEwen, J. D., Cai, X., Kitching, T. D., & LSST Dark EnergyScience Collaboration. 2019, MNRAS, 489, 3236Price, M. A., McEwen, J. D., Pratley, L., & Kitching, T. D. 2020b, arXiv e-prints,arXiv:2004.07855Repetti, A., Pereyra, M., & Wiaux, Y. 2019, SIAM Journal on Imaging Sciences,12, 87Schrabback, T., Hartlap, J., Joachimi, B., et al. 2010, A&A, 516, A63Scoville, N., Abraham, R. G., Aussel, H., et al. 2007, ApJS, 172, 38Shirasaki, M., Yoshida, N., Hamana, T., & Nishimichi, T. 2012, ApJ, 760, 45Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv preprint arXiv:1503.03757Starck, J.-L., Candès, E., & Donoho, D. 2003, A&A, 398, 785–800Starck, J.-L., Elad, M., & Donoho, D. 2004, Advances in Imaging and ElectronPhysics, 132Starck, J.-L., Moudden, Y., Abrial, P., & Nguyen, M. 2006, A&A, 446, 1191–1204Starck, J.-L., Murtagh, F., & Fadili, J. 2015, Sparse Image and Signal Processing:Wavelets and Related Geometric Multiscale Analysis (Cambridge UniversityPress)Starck, J.-L., Pires, S., & Réfrégier, A. 2006, A&A, 451 [ astro-ph/0503373 ],1139-1150Tallada, P., Carretero, J., Casals, J., et al. 2020, Astronomy and Computing, 32,100391Teyssier, R. 2002, A&A, 385 [ astro-ph/0111367 ], 337-364Upham, R. E., Whittaker, L., & Brown, M. L. 2020, MNRAS, 491, 3165Wagner-Carena, S., Hopkins, M., Diaz Rivero, A., & Dvorkin, C. 2020, MN-RAS, 494, 1507Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, Phys. Rev. D, 70,083511Wiener, N. 1949, Extrapolation, interpolation, and smoothing of stationary timeseries: with engineering applications, Vol. 7 (MIT press Cambridge)Zaroubi, S., Ho ff man, Y., Fisher, K. B., & Lahav, O. 1995, Astrophysical Journal,449, 446man, Y., Fisher, K. B., & Lahav, O. 1995, Astrophysical Journal,449, 446