Weak Lensing Mass Mapping and Higher-Order Statistics for Precision Cosmology
December 2nd, 2025
Andreas Tersenov
PhD Day, December 2nd, 2025
© Greg Bacon/STScI/NASA Goddard Space Flight Center
What do we want to know about the Universe?
The Beginning
\(A_s, n_s, f_{NL}, \tau...\)
Statistics of initial conditions
The Content
\(\Omega_m, \Omega_b, m_\nu,...\)
Cosmic ingredients
The Future
\(H_0, \Omega_\Lambda, w_0, w_a...\)
Expansion & dynamics
ΛCDM model
The simplest framework that fits almost all cosmological observations with a minimal parameter set.
But it is not without its challenges.
- How is matter distributed in the Universe?
- How does structure grow over time?
- Cosmological tensions (e.g., \(H_0\)).
Probes of the Universe
| Description | Param. | Value |
|---|---|---|
| Hubble parameter | \( H_0 \) | \( 67.66 \pm 0.42\ \mathrm{km\,s^{-1}\,Mpc^{-1}} \) |
| Total matter density | \( \Omega_m \) | \( 0.3111 \pm 0.0056 \) |
| Dark matter density | \( \Omega_c h^2 \) | \( 0.11933 \pm 0.00091 \) |
| Baryon density | \( \Omega_b h^2 \) | \( 0.02242 \pm 0.00014 \) |
| Dark energy density | \( \Omega_\Lambda \) | \( 0.6889 \pm 0.0056 \) |
| Power spectrum normalisation | \( \sigma_8 \) | \( 0.8102 \pm 0.0060 \) |
| Spectral index | \( n_s \) | \( 0.9665 \pm 0.0038 \) |
| Reionisation optical depth | \( \tau \) | \( 0.0561 \pm 0.0071 \) |
| Sum of neutrino masses | \( M_\nu \) | \( < 0.12\ \mathrm{eV} \) |
Planck Collaboration, “Planck 2018 results. VI. Cosmological parameters”, A$\&$A
Different cosmological probes capture different physical aspects of the Universe:
Clockwise: CMB - initial conditions | BAO & SNe - geometry |
Galaxy clustering - matter growth | Weak lensing - total matter
© ESA / Euclid Consortium
Convergence \( \kappa \)
As a Line-of-sight (LOS) projection of the 3D overdensity:
\[ \kappa(\boldsymbol{\theta}) = \int_{0}^{\chi_s}\! d\chi\; W(\chi)\,\delta(\chi\,\boldsymbol{\theta},\chi),\qquad W(\chi)=\frac{3H_0^2\Omega_m}{2c^2}\,\frac{\chi(\chi_s-\chi)}{a(\chi)\,\chi_s}. \]
- Encodes both geometry (via \( \chi, a(\chi) \)) and growth (via \( \delta \)).
Simulated \( \kappa(\boldsymbol{\theta}) \) map.
The Bayesian Framework for Cosmological Inference
Parameters
\( \boldsymbol{\theta} \)
\( \boldsymbol{\theta} \)
→
Likelihood
\( p(x\,|\,\boldsymbol{\theta},\mathcal{M}) \)
Forward model of data given parameters
\( p(x\,|\,\boldsymbol{\theta},\mathcal{M}) \)
Forward model of data given parameters
→
Posterior
\( p(\boldsymbol{\theta}\,|\,x,\mathcal{M}) \)
\( p(\boldsymbol{\theta}\,|\,x,\mathcal{M}) \)
\[
p(\boldsymbol{\theta}\,|\,x,\mathcal{M})
\propto
\underbrace{p(x\,|\,\boldsymbol{\theta},\mathcal{M})}_{\text{likelihood}}
\underbrace{p(\boldsymbol{\theta}\,|\,\mathcal{M})}_{\text{prior}}
\]
Interpretation
- Prior - belief about parameters before seeing data.
- Likelihood - probability of data given parameters.
- Posterior - updated belief about parameters given data.
Inference pipeline
Observed data
\( \mathbf{d}_{\mathrm{obs}} \)
\( \mathbf{d}_{\mathrm{obs}} \)
↓
Simulations / Theory
\( \mathbf{d}_{\mathrm{pred}}(\boldsymbol{\theta}) \)
\( \mathbf{d}_{\mathrm{pred}}(\boldsymbol{\theta}) \)
↓
Summary statistic
\( s(\cdot) \)
↓
Likelihood
\( p\!\left(s_{\mathrm{obs}}\mid\boldsymbol{\theta}\right) \)
\[ \mathrm{L}(\boldsymbol{\theta}) \propto \exp\!\Big[-\tfrac{1}{2}\,\Delta s^{\!\top}\,\mathbf{C}^{-1}\,\Delta s\Big], \quad \Delta s \!=\! s_{\mathrm{obs}}-s_{\mathrm{theory}}(\boldsymbol{\theta}) \]
\( p\!\left(s_{\mathrm{obs}}\mid\boldsymbol{\theta}\right) \)
\[ \mathrm{L}(\boldsymbol{\theta}) \propto \exp\!\Big[-\tfrac{1}{2}\,\Delta s^{\!\top}\,\mathbf{C}^{-1}\,\Delta s\Big], \quad \Delta s \!=\! s_{\mathrm{obs}}-s_{\mathrm{theory}}(\boldsymbol{\theta}) \]
MCMC
↓
Inference
\( p\!\left(\boldsymbol{\theta}\mid s_{\mathrm{obs}}\right) \)
\( p\!\left(\boldsymbol{\theta}\mid s_{\mathrm{obs}}\right) \)
Traditional Approach: The 2-Point Route
Concept
- Measure galaxy ellipticities \( \epsilon \) .
- Compute the 2-pt functions.
- Run MCMC with an analytic likelihood to recover posteriors:
\[ p(\theta \mid x)\;\propto\; \underbrace{p(x \mid \theta)}_{\text{likelihood}}\, \underbrace{p(\theta)}_{\text{prior}} \]
Traditional Approach: The 2-Point Route
KiDS-1000 cosmic shear constraints using 2-point functions only
(Asgari et al. 2022)
Concept
- Measure galaxy ellipticities \( \epsilon \) .
- Compute the 2-pt functions.
- Run MCMC with an analytic likelihood to recover posteriors:
\[ p(\theta \mid x)\;\propto\; \underbrace{p(x \mid \theta)}_{\text{likelihood}}\, \underbrace{p(\theta)}_{\text{prior}} \]
…but is this sufficient?
Traditional Approach: The 2-Point Route
Traditional Approach: The 2-Point Route
- Efficient, but lose non-Gaussian info.
- Non-Gaussian structure lives in higher orders.
- To access it → go beyond 2-pt: Higher-Order Statistics (HOS) (bispectrum, peak counts, Minkowski Functionals etc).
Indeed, HOS can tighten constraints significantly
Simulation-Based Inference
Le problème en question
- Weak lensing inference with higher-order statistics (HOS) often relies on mass maps.
- Producing these maps is an ill-posed inverse problem — there is no unique solution.
- Current surveys use simplistic algorithms because they are computationally efficient.
This raises an important question:
Does the choice of mass mapping algorithm affect the inferred cosmological parameters?
Or does it not matter, as long as the same operator is applied to both data and simulations?
How to test this properly?
We need to create a pipeline that can:
$\,$Introduction - Weak Lensing$\,$
- WL = Observational technique in cosmology for studying the matter distribution in the universe
- Principle: deflection of light from distant galaxies by gravitational fields → causes image distortion
- Weak → subtle & coherent distortions of background galaxy shapes
- WL provides a direct measurement of the gravitational distortion.
- WL enables us to probe the cosmic structure, investigate the nature of dark matter, and constrain cosmological parameters.
Shear & Convergence
Convergence $\kappa$
$$ \kappa = \frac{1}{2} (\partial_1 \partial_1 + \partial_2 \partial_2) \psi =
\frac12
\nabla^2 \psi $$
⟶ difficult to measure
Shear $\gamma$
$$ \gamma_1 = \frac{1}{2} (\partial_1 \partial_1 - \partial_2 \partial_2) \psi,
\,\,
\gamma_2 = \partial_1 \partial_2 \psi $$
⟶ can be measured by statistical analysis of galaxy shapes
Shear estimation
Galaxy shapes as estimators for gravitational shear
$$ e = \gamma + e_i \qquad \mbox{ with } \qquad e_i \sim \mathcal{N}(0, I)$$
- We are trying the measure the ellipticity $e$ of galaxies as an estimator for the gravitational shear $\gamma$
Relation between $\kappa$ and $\gamma$
- From convergence to shear: $\,\, \gamma_i = \hat{P}_i \kappa$
- From shear to convergence: $\,\, \kappa = \hat{P}_1 \gamma_1 + \hat{P}_2 \gamma_2$
$$ \hat{P}_1(k)=\dfrac{k_1^2 - k_2^2}{k^2}, \,\,
\hat{P}_2(k)=\dfrac{2k_1k_2}{k^2} $$
The Weak Lensing Mass-Mapping as an Inverse Problem
Shear $\gamma$
![]()
Convergence $\kappa$
![]()
$$\gamma_1 = \frac{1}{2} (\partial_1^2 - \partial_2^2) \ \Psi \quad;\quad \gamma_2 = \partial_1
\partial_2 \ \Psi \quad;\quad \kappa = \frac{1}{2} (\partial_1^2 + \partial_2^2) \ \Psi$$
$$\boxed{\gamma = \mathbf{P} \kappa}$$
Kaiser-Squires inversion
- Mass mapping problem recover an estimate of the convergence from measurements of the shear $$\kappa = \dfrac12 \Delta \psi, \,\, \gamma_1 = \dfrac12 \left( \partial_1^2 \psi - \partial_2^2 \psi \right), \,\, \gamma_2 = \partial_1 \partial_2 \psi $$ $$\Rightarrow \kappa = \Delta^{-1} \left[ \left( \partial_1^2 - \partial_2^2 \right) \gamma_1 + 2\partial_1 \partial_2 \gamma_2 \right] $$
Kaiser-Squires estimator:
$$ \tilde{\kappa}_E+i \tilde{\kappa}_B=\left(\frac{k_1^2-k_2^2}{k^2}+i \frac{2 k_1
k_2}{k^2}\right)\left(\tilde{\gamma}_1+i
\tilde{\gamma}_2\right)=\mathbf{P}\left(\tilde{\gamma}_1+i \tilde{\gamma}_2\right)
$$
→ $\tilde{\kappa}_E$ and $\tilde{\kappa}_B$ are complex representations of convergence E and B modes.
From galaxies to mass maps
Mass mapping as an inverse problem
- Binned data: $\gamma = F^* P F \, \kappa$
- Unbinned data: $\gamma = T^* P F \, \kappa$
$T =$ Non Equispaced Disctrete Fourier Transform (NDFT)
$$ \tilde{\boldsymbol{\kappa}}=\min _{\boldsymbol{\kappa}}\|\boldsymbol{\gamma}-\mathbf{A}
\boldsymbol{\kappa}\|_{\boldsymbol{\Sigma}_n}^2 $$
with $\mathbf{A} = \boldsymbol{F} \mathbf{P} \boldsymbol{F}^* $
There is no unique and stable solution, it is an ill-posed problem.
Kaiser-Squires inversion
Advantages:
- Simple linear operator
- Very easy to implement in Fourier space
- Optimal, in theory
Practical difficulties:
- Shear measurements are discrete, noisy, and irregularly sampled
- We actually measure the reduced shear: $g = \gamma/(1-\kappa)$
- Masks and integration over a subset of $\mathbb{R}^2$ lead to border errors ⇒ missing data problem
- Convergence is recoverable up to a constant ⇒ mass-sheet degeneracy problem
Inpainting
- Inpainting methods → techniques to fill in missing data in images by estimating the missing values from the known ones.
- A way to deal with the missing data problem in mass mapping and with the border effects in the convergence field.
- Implement a sparse inpainting based on the DCT (assume that the full $\kappa$-field is sparse in the DCT dictionary → apply an iterative thresholding algorithm to recover the missing data).
Bayesian reconstruction
- Mass mapping problem → statistical inference problem
- Goal: infer most probable value of $\kappa$-field given observed shear data
- Maximum A Posteriori solution: $\hat{x} = \arg \max_x \, \log p(y|x) + \log p(x)$
-
In practice: this is usually solved iteratively by
alternating two steps:
- Move toward better data fit (gradient of likelihood)
- Enforce the prior using a proximal operator
Proximal operator
- Acts as a "smart denoiser" by finding the closest solution that satisfies the prior.
- Example: sparsity prior → proximal operator performs thresholding to enforce sparsity in the solution.
Wiener filter
- Assumes prior on $\kappa$ → Gaussian random field $p_{\rm Gauss}(\kappa) = \frac{1}{\sqrt{\det 2\pi \mathbf{S}}} \exp{\left( - \frac12 \tilde{\kappa}^{\dagger} \mathbf{S}^{-1} \tilde{\kappa} \right)}$
- Likelihood (assuming uncorrelated, Gaussian noise) → also Gaussian $p(\mathbf{\gamma} | \kappa) = \frac{1}{\sqrt{2\pi\det \mathbf{N}}} \exp{\left[ -\frac12 (\mathbf{\gamma - \mathbf{A} \mathbf{\kappa}})^{\dagger} \mathbf{N}^{-1} ( \gamma - \mathbf{A} \mathbf{\kappa} ) \right]}$
- The Gaussian prior encodes the assumption that the fluctuations in the $\kappa$-field are well described by a Gaussian random field, with power spectrum given by the cosmological model
Convergence power spectrum
$$ \langle \tilde{\kappa}(\boldsymbol{k}) \tilde{\kappa}^*(\boldsymbol{k}') \rangle =
(2\pi)^2
\delta_D(\boldsymbol{k} - \boldsymbol{k}') P_{\kappa}(k) $$
Stat. measure of the spatial distribution of the convergence field → quantifies the
amplitude of the fluctuations in κ as function of their spatial scale
Wiener filter
Wiener solution of the inverse problem
$$\hat{\kappa}_{\text {wiener }}=\arg \min _\kappa\left\|\Sigma^{-1 /
2}\left(\gamma-\mathbf{F}^* \mathbf{P F} \kappa\right)\right\|_2^2+\log p_{\text {Gaussian
}}(\kappa)$$
- This solution corresponds to the maximum a posteriori (MAP) solution under the assumption of a Gaussian prior on $\kappa$, and it matches the mean of the Gaussian posterior.
Wiener reconstruction
$$\hat{\kappa}_{\text {wiener }}=\mathbf{S} \mathbf{P}^{\dagger}
\left[\mathbf{P} \mathbf{S} \mathbf{P}^{\dagger} + \mathbf{N} \right]^{-1} \tilde{\gamma}$$
Sparse recovery
- Decomposes the signal into another domain (dictionary), where it is sparse
- Implement the wavelet transform → decomposes the signal into wavelet functions (waveforms of limited duration with an average value of zero)
- Use starlet wavelets → represent well structures resembling the $\kappa$ of a DM halo (positive & isotropic)
- The application of sparsity prior enforces a cosmological model where the matter field is a combination of spherically symmetric DM halos
MCALens
- Models $\kappa$-field as a sum of a Gaussian and non-Gaussian component $$\boxed{\kappa = \underbrace{\kappa_{\rm NG}}_{\text{Standard Wiener filter approach}} + \underbrace{\kappa_\rm G}_{\text{Modified wavelet approach}}}$$ $$\min _{\kappa_G, \kappa_{N G}}\left\|\gamma-\mathbf{A}\left(\kappa_G+\kappa_{N G}\right)\right\|_{\Sigma_n}^2+C_{\mathrm{G}}\left(\kappa_G\right)+C_{\mathrm{NG}}\left(\kappa_{N G}\right) $$
- MCA (morphological Component Analysis) performs an alternating
minimization
scheme:
- Estimate $\mathbf{\kappa}_{\mathrm{G}}$ assuming $\mathbf{\kappa}_{\mathrm{NG}}$ is known: $$\min _{\kappa_G}\left\|\left(\gamma-\mathbf{A}\kappa_{N G}\right) -\mathbf{A}\kappa_{G} \right\|_{\Sigma_n}^2+C_{\mathrm{G}}\left(\kappa_G\right)$$
- Estimate $\mathbf{\kappa}_{\mathrm{NG}}$ assuming $\mathbf{\kappa}_{\mathrm{G}}$ is known: $$\min _{\kappa_{NG}}\left\|\left(\gamma-\mathbf{A}\kappa_{G}\right) -\mathbf{A}\kappa_{NG} \right\|_{\Sigma_n}^2+C_{\mathrm{NG}}\left(\kappa_{NG}\right)$$
Overview of mass mapping methods
- Wiener filter: Prior on $\kappa$ → Gaussian random field, Likelihood → Gaussian, Solution → corresponds to MAP for gaussian $\kappa$
- Sparse recovery: Prior on $\kappa$ → sparse in wavelet basis, Enforces model where matter field is a combination of DM halos
- MCALens: Models $\kappa$-field as sum of a Gaussian and a non-Gaussian component, Uses Morphological Component Analysis
- Deep Learning Methods:
- DeepMass: CNN with a U-Net-based architecture, prior from simulations
- DeepPosterior: Probabilistic mass mapping with deep generative models, Prior from 2pt statistics modelling at large scales & Deep Learning on simulations for small scales, Sampling with Annealed HMC
Does the choice of mass-mapping method matter for cosmology?
A. Tersenov, L. Baumont, J.L. Starck, M. Kilbinger, doi.org/10.1051/0004-6361/202553707
Mass mapping methods:
\[
\begin{array}{ll}
\hline \hline
\text{Method} & \text{RMSE} \downarrow \\
\hline
\text{KS } & 1.1 \times 10^{-2} \\
\text{iKS } & 1.1 \times 10^{-2} \\
\text{MCALens} & 9.8 \times 10^{-3} \\
\hline
\end{array}
\]
Higher Order Statistics: Peak Counts
- Peaks: local maxima of the SNR field $\nu = \dfrac{\left(\mathcal{W} \ast \kappa \right)(\theta_{\rm ker})}{\sigma_n^{\rm filt}}$
- Peaks trace regions where the value of $\kappa$ is high → they are associated to massive structures
Peak funtions
- The peak function is the number of peaks in a map as a function of the SNR $\nu$
Wavelet peaks
- We consider a multi-scale analysis compared to a single-scale analysis
- Apply (instead of Gaussian filter) a starlet transform → allows us to represent an image $I$ as a sum of wavelet coefficient images and a coarse resolution image
- Allows for the simultaneous processing of data at different scales → efficiency
- Each wavelet band covers a different frequency range, which leads to an almost diagonal peak count covariance matrix
Wavelet peak funtions
Starlet $\ell_1$-norm
- New higher order summary statistic for weak lensing observables
- Provides a fast multi-scale calculation of the full void and peak distribution
- $\ell_1$-norm = sum of absolute values of the starlet decomposition coefficients of a weak lensing map
- Information encoded in all pixels
- Automatically includes peaks and voids
- Multi-scale approach
- Avoids the problem of defining peaks and voids
Inference with HOS
- Use $\texttt{cosmoSLICS}$ simulations: suite designed for the analysis of WL data beyond the standard 2pt cosmic shear
- $\texttt{cosmoSLICS}$ cover a wide parameter space in $\left[ \Omega_m, \sigma_8, w_0, h \right]$.
- For Bayesian inference → use a Gaussian likelihood for a cosmology independent covariance, and a flat prior.
- To have a prediction of each HOS given a new set of parameters→ employ an interpolation with Gaussian Process Regressor (GPR)
Noise & Covariance
- We consider Gaussian, but non-white noise → the noise depends on the number of galaxies in each pixel $$\sigma_n = \dfrac{\sigma_{\epsilon}}{\sqrt{2 n_{\rm gal} A_{\rm pix}}}$$
- We incorporate masks
- Calculate covariance: $$ \boxed{C_{ij} = \sum_{r=1}^N \dfrac{\left( x_i^r - \mu_i \right) \left( x_j^r - \mu_j \right)}{N-1}} \, , \,\,\,\mu_i = \dfrac{1}{N} \sum_r x_i^r $$
From data to contours
- $\mu$: expected theoretical prediction, $d$: data array (mean over realizations of a HOS), $C$: covariance matrix
So does the choice of the mass mapping algorithm matter for the final constraints?
The (standard) mono-scale peak counts
\[
\begin{array}{lcccc}
\hline \hline
\text{FoM} & \text{KS} & \text{iKS} & \text{MCALens} \\
\hline \hline
\left(\Omega_m, h\right) & 476 & 453 & 450 \\
\left(\Omega_m, w_0\right) & 152 & 141 & 233 \\
\left(\Omega_m, \sigma_8\right) & 1323 & 1285 & 1740 \\
\left(h, w_0\right) & 55 & 63 & 87 \\
\left(h, \sigma_8\right) & 336 & 292 & 293 \\
\left(w_0, \sigma_8\right) & 75 & 72 & 124 \\
\hline
\left(\Omega_m, h, w_0, \sigma_8\right) & 492 & 444 & 578 \\
\hline
\end{array}
\]
Wavelet multi-scale peak counts
\[
\begin{array}{lcccc}
\hline \hline
\text{FoM} & \text{KS} & \text{iKS} & \text{MCALens} \\
\hline \hline
\left(\Omega_m, h\right) & 670 & 702 & 2159 \\
\left(\Omega_m, w_0\right) & 247 & 244 & 1051 \\
\left(\Omega_m, \sigma_8\right) & 2414 & 2517 & 9039 \\
\left(h, w_0\right) & 82 & 80 & 259 \\
\left(h, \sigma_8\right) & 411 & 433 & 1335 \\
\left(w_0, \sigma_8\right) & 131 & 129 & 577 \\
\hline
\left(\Omega_m, h, w_0, \sigma_8\right) & 758 & 755 & 1947 \\
\hline
\end{array}
\]
Where does this improvement come from?
Summary statistics
Part 2
A plug-and-play approach with fast uncertainty quantification for weak lensing mass mapping
H. Leterme, A. Tersenov, J. Fadili, and J.-L. Starck (in prep.)
Wait... what about Deep Learning for Mass Mapping?
Yep, people have tried it! ...And it works!
Example: DeepMass
So what's the problem?
Plug-and-Play Mass Mapping
Main idea:
- Use PnP framework: replace prox by an on-the-shelf deep denoiser trained on simulations $$ \kappa^{n+1} = \mathrm{prox}_{\tau g} \left( \kappa^n - \tau \nabla f(\kappa^n) \right) $$ $$ \kappa^{n+1} = F_{\Theta} \left( \kappa^n - \tau \mathbf{B} \left( \mathbf{A} \kappa^n - \gamma \right) \right) $$ $$ \kappa^{n+1} = \mathrm{Denoiser} \left( \kappa^n + \mathrm{Data \, residual} \right) $$
- Series converges towards a fixed point $\hat{\kappa}$
- If we choose $\mathbf{B} = \mathbf{A}^T \Sigma^{-1}$ → training phase independent of the noise covariance matrix
Instead of explicitly writing down a prior, we learn what a "likely" $\kappa$ looks like from simulations and enforce it through denoising.
Implementation
Training
- Denoiser models implemented: DRUNet & SUNet
- Trained on κTNG and cosmoSLICS simulations, using pairs of $(\kappa_{\rm true}, \gamma_{\rm obs})$ as training data
Uncertainty bounds
