PhD defense, University of Crete, 14 September 2026

Trustworthy non-Gaussian inference for weak-lensing cosmology

Doctoral thesis in Physics

slides atandreastersenov.github.io/talks/PhD_Defense_2026

The Universe evolved from quantum fluctuations to galaxies over 13.8 billion years

history of the Universe: inflation, the cosmic microwave background, the dark ages, the growth of structure, and accelerated expansion

ΛCDM describes the Universe with six free parameters, plus a set of assumptions

Gravity: general relativity, unmodified
Geometry: spatially flat, homogeneous and isotropic on large scales
Contents: cold dark matter (collisionless, non-relativistic), baryons, radiation, and a constant $\Lambda$
Initial conditions: adiabatic, near-Gaussian and near-scale-invariant, from inflation

energy budget of the Universe: about 70 per cent dark energy and 30 per cent total matter, with the matter slice broken down further

history of the Universe, annotated with the three things ΛCDM does not explain

Independent probes test the model through both geometry and growth

four cosmological probes: the cosmic microwave background, baryon acoustic oscillations with supernovae, galaxy clustering, and weak gravitational lensing

Foreground mass deflects the light of background galaxies and distorts their images

light from three background galaxies travelling through spacetime curved by foreground dark matter, arriving distorted in the observed image

Image © American Physical Society / Alan Stonebraker; galaxy images from STScI/AURA, NASA, ESA and the Hubble Heritage Team / astronomie.nl.

Euclid: Mapping the Dark Universe

ESA mission, launched 2023

Lensing and clustering together, to measure cosmic geometry and growth at percent-level precision.

14,000 deg², about a third of the sky
The shapes of billions of galaxies
Weak lensing from the visible imaging, clustering from the near-infrared
First cosmological release, June 2027

ESA / Euclid Consortium

Real weak-lensing cosmological analyses are extremely complicated

the DES Year 3 three-by-two-point analysis flowchart: wide-field images and deep-field photometry feeding PSF modelling, shape catalogues, image simulations, blending and redshift calibration, then two-point measurements, analysis choices and modelling, and finally the cosmology contours — DES Y3 3×2pt — figure: Alexandra Amon, DES collaboration

Almost every box is a paper on its own, and most of the work is measurement and calibration (PSF modelling, shape calibration, blending, redshift distributions, covariances). .

The analysis chain from galaxy shapes to parameters

the first half — papers 1 and 2

Making the mass map from the measured shear.

the second half — papers 3 and 4

Reading the map, from a summary statistic to the posterior.

The analysis chain from galaxy shapes to parameters, and the four questions of this thesis

measured galaxy shapes the data

→

reconstruct the matter map step 1

→

summarise it with a statistic step 2

→

posterior contours on the cosmological parameters

infer the parameters step 3

Every one of those steps is now built out of learned components, and every one of them can bias the result or discard information, and there is no internal check that would detect it.

The next two parts live in one step of that chain

Shear in, a mass map out. Everything in Parts 1 and 2 is about what happens between those two boxes.

Part 1: Impact of weak-lensing mass-mapping algorithms on cosmology inference

A. Tersenov, L. Baumont, J.L. Starck, M. Kilbinger, arXiv:2501.06961

We measure the shear, but theory predicts the convergence

three panels showing how a circular source is distorted: convergence enlarges it isotropically, and the two shear components stretch it along the axes and along the diagonals

Convergence $\kappa$

$$\kappa \;=\; \tfrac12\left(\partial_1\partial_1 + \partial_2\partial_2\right)\psi \;=\; \tfrac12\nabla^2\psi$$

isotropic magnification
proportional to the projected mass
not directly observable (the unlensed size is unknown)

Shear $\gamma$

$$\gamma_1 = \tfrac12\left(\partial_1\partial_1 - \partial_2\partial_2\right)\psi, \quad \gamma_2 = \partial_1\partial_2\psi$$

anisotropic distortion
coherent across neighbouring galaxies
measurable from galaxy ellipticities

The convergence is a projection of the matter density, sensitive to geometry and growth

Why we want the convergence

it is a scalar field, where the shear is spin‑2
it is the projected matter along the line of sight, weighted by how efficiently each piece of it lenses

So a convergence map is a map of the mass.

the convergence map of the UNIONS survey footprint, reconstructed with MCALens: filamentary structure across two disconnected sky patches, with the survey mask left grey — Hervas‑Peters, …, Tersenov et al. 2026

Shear and convergence are both second derivatives of the lensing potential

a shear map and a mass map, each linked to the lensing potential in the middle by a pair of second-derivative relations

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$
linear and exact for a complete, noiseless field

Relation between $\kappa$ and $\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$

In theory optimal.
In practice, $\{ \textrm{noise, masks, irregular sampling of galaxies, mass-sheet degeneracy} \} \longrightarrow$ ill-posed inverse problem

In practice, mass mapping is an ill-posed inverse problem

$$\gamma \;=\; \underbrace{\mathbf{M}}_{\text{the mask}}\quad \underbrace{\mathbf{P}}_{\text{shear operator}}\; \kappa \;+\; \underbrace{\mathbf{n}}_{\text{shape noise}}$$

$\kappa = \;?$

a smooth, heavily denoised convergence map

one of many convergence maps consistent with the measurement

all consistent with the data

shape noise $\mathbf{n}$, intrinsic ellipticities much larger than the shear
survey mask $\mathbf{M}$, missing data and a non-local inversion
mass-sheet degeneracy, a constant added to $\kappa$ leaves the shear unchanged

Selecting one solution requires a prior on $\kappa$, and the choice of prior is what distinguishes the methods.

Kaiser–Squires inversion is exact for complete, noiseless data, but amplifies the noise

the true convergence map beside its Kaiser-Squires reconstruction, the structure of the first buried in noise in the second

Advantages

linear, no free parameters
one FFT, fast on any survey
exact on a complete, noiseless field

Limitations

amplifies shape noise at small scales
mask effects leak across the map (non-local operator)

Every mass-mapping method minimises a data term plus a regulariser that encodes the prior

Find the $\kappa$ that fits the measured shear and is plausible under a prior

$$\hat{\kappa} \;\in\; \arg\min_{\kappa}\; \underbrace{\tfrac{1}{2}\big\lVert \mathbf{P}\kappa - \gamma \big\rVert^{2}_{\mathbf{N}^{-1}}}_{\text{data fidelity: does it match the shear?}} \;+\; \lambda\, \underbrace{\mathcal{R}(\kappa)}_{\text{regulariser: is it a plausible map?}}$$

$\mathbf{N}$, noise covariance
$\lambda$, regularisation strength
the data term is the same for all methods; they differ in $\mathcal{R}$

MCALens models κ as a Gaussian plus a sparse component

$\kappa \;=$

$\kappa_{\rm G}$ Gaussian component

estimated with a Wiener filter, which needs only the power spectrum.

$+$

$\kappa_{\rm NG}$ non-Gaussian component

sparse in the starlet domain (peaks and haloes).

Alternating minimisation

solve $\kappa_{\rm G}$, holding $\kappa_{\rm NG}$ fixed ↻ solve $\kappa_{\rm NG}$, holding $\kappa_{\rm G}$ fixed

Iterate to convergence. Each sub-problem is solved with a proximal step.

the proximal operator

$$\mathrm{prox}_{\lambda\mathcal{R}}(v) \;=\; \arg\min_{u}\; \tfrac{1}{2}\lVert u - v \rVert^{2} \;+\; \lambda\,\mathcal{R}(u)$$

a gradient step on the data term gives $v$
prox returns the point closest to $v$ that is consistent with the prior

Different priors give different maps. Does the difference affect the cosmological results?

A practical question for Euclid: is an advanced reconstruction method worth the effort, or is any reasonable method good enough?

§1 The simulations, statistic and likelihood are fixed, and only the mass-mapping method varies

simulations

cosmo‑SLICS

shear maps

with Euclid‑like observational effects

mass maps

KS, iKS, MCALens

the only thing that varies

summary statistic

set of non‑Gaussian statistics

$\log \mathcal{L}(\theta)$

inference

posterior

Ω_m, σ₈, h, w₀

the output we compare

Everything else is fixed, so any change in the posterior comes from the reconstruction.

§1 MCALens recovers small-scale structure that Kaiser–Squires smooths out

four panels of the same simulated convergence field: the truth, Kaiser-Squires, iterative Kaiser-Squires and MCALens, with the survey mask outlined

prior on $\kappa$

truth the true map (simulations only)

KSEuclid baseline none. Direct inversion, then smoothing

iKSEuclid baseline none; the mask is inpainted (DCT) before inversion

MCALensstate of the art a Gaussian background plus a sparse non-Gaussian component

Mass mapping has been treated as a preprocessing step with little effect on the cosmology, which is why the Euclid baseline is the simplest method, Kaiser–Squires.

§1 MCALens improves the figure of merit by 157 % over Kaiser–Squires, by making more information in the data available.

four-parameter corner plot comparing the posteriors from KS, iKS and MCALens wavelet peak counts

the four-parameter figure of merit relative to Kaiser-Squires, drawn as stems: Kaiser-Squires and inpainting KS both at 1.00, MCALens at 2.6; the colours are the corner plot's

MCALens improves the RMSE by only 4 %, but the figure of merit by 157 %
⇒ WL surveys should use advanced reconstruction methods to maximise their scientific outputs.

§1 The gain comes from scales below 8′, which only MCALens recovers

corner plot for Kaiser-Squires with three, four, five and six starlet scales included; the contours stop shrinking after the fourth — adding the 4′ and 2′ scales does not change the contours

the same corner plot for MCALens; the contours continue to shrink as each finer scale is added — every added scale tightens the contours

interpretation MCALens recovers small‑scale structure that the linear inversion loses to noise.

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
arXiv:2603.22006

What about Deep Learning for Mass Mapping?

Yes, people have tried it (with various approaches) … and it works!

However…

mass mapping method	type	accurate	flexible	fast rec.	fast UQ
Iterative Wiener	model-driven, Gaussian prior	✗	✓	✓	✗
MCALens	model-driven, Gaussian + sparse	≈	✓	✗	✗
DeepMass	data-driven, U-Net	✓	✗^*	✓	✓
DeepPosterior	data-driven, U-Net + MCMC	✓	✓	✗	✗
MMGAN	data-driven, GAN	✓	✗^*	≈	≈
what we’d like	data-driven	✓	✓	✓	✓

^* must be retrained when the noise level or the mask changes

Deep Learning for Mass Mapping?

Mass mapping method	Type	Accurate	Flexible	Fast rec.	Fast UQ
Iterative Wiener	Model-driven (Gaus. prior)	✗	✓	✓	✗
MCALens	Model-driven (Gaus. + sparse)	≈	✓	✗	✗
DeepMass	Data-driven (UNet)	✓	✗*	✓	✓
DeepPosterior	Data-driven (UNet + MCMC)	✓	✓	✗	✗
MMGAN	Data-driven (GAN)	✓	✗*	≈	≈
What we'd like	Data-driven	✓	✓	✓	✓

Plug-and-Play Mass Mapping

PnPMass, iteration 1: the measured shear, initialisation at zero, the forward step giving a noisy map, and the denoiser giving a clean one beside the ground truth

PnPMass, iteration 2: the output of the first iteration fed back in

Use the PnP framework: replace prox by an off‑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)$

The series converges to a fixed point $\hat{\kappa}$
choosing $\mathbf{B} = \mathbf{A}^T \Sigma^{-1}$ makes the method indep. of the noise covariance

plug-and-play

from the optimisation literature; inverse problem with a learned prior, very flexible
start from the Part 1 iteration: gradient step on the data term, then proximal step
PnP replaces the proximal operator with a deep denoiser trained on simulations → converges to a fixed point
better than analytical priors: the denoiser learns a far more complex, realistic prior
Swin transformer, 7 million parameters, trained once on white Gaussian noise over a range of levels
▲ noise covariance and mask enter only in the data term, at inference → same network for any noise, any footprint. That is the flexibility
walk-through: two shear components, κ initialised at zero → forward (gradient) gives a noisy map → backward (denoiser) projects onto realistic mass maps
» output fed back, iterate to convergence
short path if behind: eight iterations, alternating data step and denoiser; residual variant slightly better, in the backup

How do we estimate uncertainties?

Step 1

We train a second neural network to estimate the posterior variance of $\kappa$: $$ \arg\min_{\Omega} \mathbb{E} \left[ \left\| G_{\Omega}(\gamma) - \big(\kappa - F_{\Theta}(\gamma) \big)^2 \right\|^2_2 \right] $$
Trained on simulated pairs $(\kappa, \gamma)$, just like the denoiser — but now focused on uncertainty.
This gives fast, pixel-wise error bars
Uncertainty estimation is fast — just one extra iteration after κ̂ is computed → adds almost no overhead

Step 2

Neural networks tend to produce miscalibrated uncertainties.
We apply conformal quantile regression (CQR) to adjust uncertainty intervals so they have guaranteed statistical coverage.
CQR uses a held-out calibration set to compute a correction factor for each pixel.
Result: Reliable, data-driven uncertainty maps — with built-in coverage guarantees.

§2 PnPMass is as accurate as the fine‑tuned networks, with the smallest calibrated error bars

mean calibrated interval length against reconstruction error for Wiener, MCALens, DeepMass and both PnPMass variants — 512 test maps, κTNG. Circles before calibration, diamonds after.

lower left is better; after calibration, all methods reach the target coverage
accuracy comparable to the SOTA deep‑learning methods
smallest calibrated intervals of all methods
PnPMass is indep. of mask and noise covariance

the uncertainty map traces the structure

a PnPMass convergence map beside its per-pixel calibrated uncertainty, on the same field — The map, and its calibrated σ on the same field.

DeepMass and MMGAN need retraining for every new mask or noise level; PnPMass is trained once and applies to any configuration.

Performance

Work in progress: extending to tomographic mass mapping, and full 3D mass reconstruction.

Conclusions

Investigated how different mass mapping algorithms affect cosmological inference using HOS from WL data.
Constructed new mass mapping algorithm based on the PnP formalism.

Results

With a state-of-the-art mass-mapping method (MCALens) we managed to get $\sim 157\%$ improvement in FoM over KS.
Increase in constraining power comes from the more accurate recovery of the smaller scales.
Wavelet Peak Counts: Provide tighter constraints than single-scale peak counts.
PnP Mass Mapping:

Provides a fast, flexible, and accurate mass mapping algorithm that can be used with any denoiser.
Provides fast and reliable uncertainties using moment networks and conformal quantile regression.

Takeaway

Mass-mapping Matters: Choosing an advanced mass mapping method significantly enhances constraints on cosmological parameters from HOS.

superseded — not part of the talk

LAM originals, Part 1

Replaced on 2026-09-03 by the restyled Act 1 slides. Kept for reference; nothing is deleted.

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} \]

Inference pipeline

Observed data
$ \mathbf{d}_{\mathrm{obs}} $

↓

Simulations / Theory
$ \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}) \]

MCMC

↓

Inference
$ p\!\left(\boldsymbol{\theta}\mid s_{\mathrm{obs}}\right) $

So does the choice of the mass mapping algorithm matter for the final constraints?

Inference results

Mono-scale peak counts

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?

alternatives — not part of the talk

Restyled Acts 1–2

The same argument rebuilt in the preprint components. Kept for comparison against the lifted originals.

act 1 — thesis chapter 2

Is mass mapping preprocessing, or does the choice of reconstruction change the cosmology we infer?

Tersenov, Baumont, Starck & Kilbinger, arXiv:2501.06961

QR code linking to the paper on the impact of mass mapping on cosmological inference

act 2 — thesis chapter 3

Can a reconstruction be flexible, fast, accurate and honest about its own uncertainty at the same time?

Leterme, Tersenov, Starck et al. — joint work: I co‑developed the method and built the uncertainty quantification

§2 Every existing reconstruction gives up at least one of the four

model‑driven

Wiener, MCALens

flexible — the physics is written down, so nothing is retrained

no useful uncertainty; Wiener assumes the field is Gaussian; MCALens is slow

data‑driven

DeepMass, MMGAN

accurate, and fast once trained

trained for one noise level and one mask; returns a point estimate with no error bar

what we want

all four at once

flexible fast accurate and calibrated

—

§2 PnPMass: one denoiser, trained once, inside a fixed‑point iteration

the plug-and-play forward-backward loop: a gradient step towards consistency with the measured shear, followed by a learned denoising step, iterated to a fixed point — forward–backward splitting, with the hand‑designed regulariser replaced by a **learned denoiser**

$$\kappa^{\,n+1} \;=\; \mathrm{Denoiser}\big(\kappa^{\,n} + \text{data residual}\big)$$

why this buys flexibility The denoiser never sees the mask, never sees the survey, never sees a map with one particular noise level baked in. It learns one thing — what a plausible convergence field looks like. The physics of the observation stays in the gradient step, where it belongs.

Swin‑Transformer denoiser, 7.2 M parameters, trained once on white Gaussian noise with the noise level as an input channel. Eight iterations per map.

§2 As accurate as a network retrained for the observation — with the smallest calibrated error bars of any method tested

reconstructed convergence map with its per-pixel lower and upper conformal bounds — the reconstruction with its **per-pixel conformal bounds** — a lower and an upper map, not a single number

miscoverage rate per method before and after conformal calibration; afterwards every method sits on the target rate — after calibration every method hits the **same** target rate

normalised error-bar size against reconstruction RMSE for each method; PnPMass sits lowest in error-bar size — …and PnPMass gets there with the **smallest bars**

accuracy Within 1 % of DeepMass in RMSE — and DeepMass was retrained for this mask and this noise. The residual variant: within 0.5 %.

the axis that matters Conformal calibration makes every method reach the target coverage. So the comparison is which method gets there with the tightest intervals — PnPMass, on all 512 test images.

Coverage is marginal, not conditional — the miscoverage concentrates at the peaks, exactly where Chapter 2 says the information lives. And this is a single cosmology.

second half — the summaries

Parts 3 and 4 move from the map to the summary statistic

Parts 1 and 2 were about the mass map. Parts 3 and 4 are about the summary statistic, and its robustness to systematics.

The map has to be compressed into a summary statistic before inference

The map has to be compressed into a summary statistic before inference, and the choice of statistic determines how much information is retained.

The two‑point function is the standard statistic, and is complete for a Gaussian field

DES Y3 cosmic shear constraints in the Omega_m, sigma_8, S_8 plane — DES Y3 cosmic shear, two‑point functions only: harmonic ($C_\ell$) and real space ($\xi_\pm$)

what it measures

How correlated the shear is between pairs of galaxies separated by an angle $\theta$. Measured as $\xi_\pm(\theta)$, or as $C_\ell$.

every wide lensing survey measures cosmic shear this way: KiDS, DES, HSC, UNIONS
predicted analytically; its covariance is understood
two decades of work on its systematics

For a Gaussian random field, the power spectrum contains all the information. But the late‑time density field is not Gaussian.

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?

Two very different fields can have the same power spectrum

a large-scale-structure map, a Gaussian map, and their identical power spectra

To access the information the power spectrum misses we need higher‑order statistics: peak counts, wavelet statistics, the ℓ₁‑norm, Minkowski functionals.

Forecasts show that higher-order statistics constrain better than the power spectrum

forecast contours on the summed neutrino mass, the matter density and the primordial amplitude, for four summaries: the power spectrum is much the widest, peaks after a single Gaussian filter are tighter, and starlet peaks and multi-Gaussian peaks are tighter again — Forecasts from the same convergence maps: power spectrum, single-scale peaks, and multi-scale peaks. Figure: Ajani, Peel, Pettorino, Starck, Li & Liu (2021), arXiv:2001.10993

Indeed, HOS can tighten constraints significantly

Peak funtions

The peak function is the number of peaks in a map as a function of the SNR $\nu$

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

Wavelet peaks

We consider a multi-scale analysis compared to a single-scale analysis
Apply 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

$$ \boxed{ \ell_1^{j,i} = \sum_{u=1}^{\textrm{\# coef}(\mathcal{S}_{i,j})} \left| \mathcal{S}_{j,i} [u] \right| = \| \mathcal{S}_{j,i} \|_1} \, , \,\, \mathcal{S}_{j,i} = \{ w_{j,k}/B_i w_{j,k} B_{i+1} \}$$

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

Summary statistics

Key concept: wavelets

a wavelet $\psi$ is a localised, oscillating function with zero mean, $\int \psi\,\mathrm{d}x = 0$
dilations and translations of $\psi$ form a family, $\psi_{a,b}(x) = \tfrac{1}{\sqrt{a}}\,\psi\!\big(\tfrac{x-b}{a}\big)$
the coefficient $W(a,b) = \int f\,\psi_{a,b}\,\mathrm{d}x$ measures the structure at scale $a$ and position $b$

the starlet wavelet at three dilations, visibly the same shape at three widths

a signal carrying three wave packets of different wavelength at different places, decomposed into four wavelet bands plus a coarse residual, each packet answering in the band that matches its size

Wavelet peak counts and the $\ell_1$-norm are one-point statistics of the starlet coefficients

the starlet transform: a convergence map written as the sum of four wavelet-coefficient bands plus a coarse map, each band carrying structure of a characteristic angular size — the **starlet transform**: one map as a sum of band-pass images, each carrying structure of a characteristic angular size

wavelet peak counts

a smoothed convergence map with its peaks circled

local maxima of the $\kappa$ field, counted per $\kappa$ bin → focus only on a subset of the field

starlet ℓ1-norm

$$\ell_1^{\,j,i} \;=\; \sum_{u=1}^{\#\mathcal{S}_{j,i}} \big| \mathcal{S}_{j,i}[u] \big| \;=\; \lVert \mathcal{S}_{j,i} \rVert_1$$

$\mathcal{S}_{j,i}$ are the coefficients of scale $j$ whose $\kappa$ falls in bin $i$

sum of absolute coefficients per band and $\kappa$ bin → encodes information from every pixel (peaks, voids, ...)

The summary statistic is compared with theoretical predictions to infer the cosmological parameters, in a Bayesian framework

Classical inference needs an explicit likelihood, usually Gaussian in the data vector

Bayes' theorem: the probability of the parameters given the data is proportional to the probability of the data given the parameters, times the prior.

after the data

$p(\theta \mid d)$

posterior

∝

how well the parameters explain the data

$p(d \mid \theta)$

likelihood

before the data

$p(\theta)$

prior

The likelihood is usually assumed Gaussian in the data vector:

$$p(d \mid \theta) \;\propto\; \exp\!\left[-\tfrac{1}{2}\big(d - \underbrace{m(\theta)}_{\text{theory}}\big)^{\!\top} \underbrace{\mathbf{C}^{-1}}_{\text{covariance}} \big(d - m(\theta)\big)\right]$$

The posterior is then sampled with MCMC.

a Metropolis-Hastings chain wandering across a correlated Gaussian target and filling in the contours

Key concept: generative modelling

Learn the distribution from which a set of samples was drawn.

three panels: the true two-moons density, a finite set of samples drawn from it, and the density a model fits to those samples

the true distribution $\mathbb{P}$ is unknown; we observe samples $X = \{x_0, x_1, \ldots, x_n\}$ and fit a parametric model $\mathbb{P}_\theta$
a trained model can generate new samples, and evaluate the density $p_\theta(x)$
the framework does not depend on what $x$ is: the same machinery generates faces, images, video, molecules,...

A normalizing flow maps a simple distribution to a complex one through an invertible transform

a chain of invertible maps carrying a simple density through intermediate stages into a complex multi-modal one

Start from a unit Gaussian and learn an invertible map to the target distribution. The density follows from the change of variables formula:

$$\log q_\phi(x) \;=\; \log p\big(f^{-1}(x)\big) \;+\; \log\left|\det \frac{\partial f^{-1}}{\partial x}\right|$$

samples being carried by the flow, rearranging from a Gaussian cloud into the target shape — samples transported by the flow

the result a flexible model that is easy to sample from and whose density we can evaluate

Simulation‑based inference learns the posterior directly from simulated data

The likelihood $p(x \mid \theta)$ cannot be written down, but the simulator samples from it. The posterior is learned directly from simulated $(\theta, x)$ pairs.

stage one: parameters drawn from the prior, pushed through a simulator, producing simulated data

stages two and three added: the parameter-data pairs training a conditional normalizing flow, and the trained flow evaluated at the real observation to give a posterior

Part 3

The joint wavelet ℓ₁‑norm matches neural compression for tomographic weak‑lensing inference

A. Tersenov, J.-L. Starck, and M. Kilbinger

How much of the map's information can a summary statistic keep?

→

summary

→

cosmology

Beating the power spectrum is "easy". The question is how close to all of it we can get.

A neural compressor trained to be information-optimal

simulator

$\theta\sim$ prior

→

map $x$

→

network

$f_\phi$

→

summary

$t$

→

flow

$q_\psi(\theta\mid t)$

→

posterior

$p(\theta\mid t)$

\[\max_{\phi,\psi}\; I(t;\theta)\;=\;\max_{\phi,\psi}\;\mathbb{E}\,\log q_\psi\!\big(\theta\mid f_\phi(x)\big)\]

parameter space $\theta$

§3 Both summaries are compared on the same maps, through the same flow, both calibrated

kappa maps to either l1-norm or CNN-VMIM, into the same flow density estimator, to a calibration-gated posterior

§3 The CNN beats the ℓ1-norm by 36% in the figure of merit

matched-pipeline posterior for the auto-only l1-norm, with its figure of merit at 2448 — same maps, same flow, both calibrated; only the summary differs

the CNN posterior added: tighter, with a figure of merit of 3326 against the l1-norm's 2448 — same maps, same flow, both calibrated; only the summary differs

the gap

The learned summary is 36% ahead in FoM₃
But the comparison is not symmetric: something important about the maps is not taken into account

§3 But could we do better than that? Weak lensing tomography

the source galaxies as tomographic convergence maps, one per redshift bin

Credit: Justine Zeghal

§3 The cross-bin information can be reached through cross maps or through a joint statistic

route 1 — change the input

the convergence map of tomographic bin 1 — $\kappa_1$

the convergence map of tomographic bin 2 — $\kappa_2$

the pointwise product of the two maps: near-zero almost everywhere, bright where both bins have structure — $\kappa_1\kappa_2$

The pixel-wise product of two bins is strong only where both have structure. One map per bin pair, each analysed with the same ℓ1-norm.

route 2 — change the statistic

bin j coefficient the joint plane of wavelet coefficients for a bin pair, with the two one-dimensional histograms along its edges

the per-bin ℓ1-norm
only sees the two
marginal histograms bin i coefficient

$$L^{ij}_{ab} \;=\!\!\sum_{\boldsymbol{p}\,\in\,\mathrm{cell}(a,b)}\!\! \tfrac{1}{2}\big(|u_i(\boldsymbol{p})| + |u_j(\boldsymbol{p})|\big)$$

each cell sums the ℓ1 weight of the pixels that fall in it

two routes

two-point extends easily — cross-bin correlations. For HOS there is no obvious equivalent
first: extra maps for the cross structure — multiply each pair of bins pix by pix; strong only where both have structure at the same place. Compute the ℓ1 on those too
second: change the stat — define stat that generalises the ℓ1 from one map to a pair of tomo bins
each pixel has a wavelet coeff in each of the 2 maps → instead of building one l1 histogram per map (the two curves on the axes) → build one two-dimensional histogram over the pair
diagonal cells = collect places where both maps equally strong; off-diagonal = one strong, one weak
Named this stat "joint l1-norm" & you can think about it as a joint multiscale PDF in wavelet domain

§3 With joint reading of the bins, the ℓ1-norm matches the optimal CNN

matched-pipeline posteriors: the auto-only l1 first, then plus product cross-maps, then the joint l1, then the CNN landing on top of it; the FoM3 bar inset fills in alongside

Part 4

Mitigating baryonic effects in weak lensing with higher‑order statistics

A. Tersenov, S. Guerrini, J.-L. Starck, and M. Kilbinger
arXiv:2609.09131

Non-Gaussian information and baryonic feedback both live at small scales

To be able to trust our HOS results we must investigate how they are affected by each systematic at the contour level.

Illustris: baryonic feedback reshaping the cosmic web

question 1

How much does unmodelled feedback bias our higher-order statistics, and how does that grow with survey area?

question 2

Once those scales are cut, is there any constraining power left to gain over the power spectrum?

systematics = phys/obs effects in the measurement, not in the model. Unmodelled → biased result
dominant astrophysical one at small scales: baryonic feedback — AGN and supernovae push gas out of haloes, suppressing small scales
our inference sims are DMO, and feedback models disagree → physics we can't model reliably
the problem: non-Gaussian information and most systematics live on the same small scales → can't model them reliably for HOS -> conservative answer is to cut them
» optimistic: only the smallest scales contaminated, most of the non-gaus info survives
» pessimistic: contamination reaches further, and after the cut the power spectrum does just as well
two questions: how much does unmodelled feedback bias HOS, and how does it grow with area as surveys expand
And after the cut is there anything left over the power spectrum, or our efforts in designing better stats were in vein, cause they are useless in real world?

In the SBI pipeline, the statistics are computed separately on each wavelet scale

wavelet-scale maps and summary statistics

+ noise
wavelet transform

condition

Gaussian$\mathcal{N}(0,\mathbf{1})$

Density estimatorconditional MAF

Posterior$p(\theta\mid x)$

training objective

$\mathcal{L}=-\log p_\phi(\theta\mid x)$

§4 Baryonic bias grows with survey area

baryon tension in sigma versus survey area for the power spectrum, peak counts and the l1-norm: all three rise with area, the two higher-order statistics well above the power spectrum

At 14,000 deg² (Stage IV): $C_\ell$ shows 2.2σ, peaks and the $\ell_1$-norm 3.6σ
At full sky: $C_\ell$ reaches $\sim3.5\sigma$, both HOS exceed 6σ
HOS show higher bias than $C_\ell$: more sensitive to the baryonically-contaminated small scales

baryon-biased power-spectrum posteriors, tightening and marching away from the truth as the survey area grows from 2,000 to 35,000 square degrees

measured at full map resolution: $\ell_{\rm max}=1024$ for $C_\ell$, all four wavelet bands for the HOS

§4 Removing the bias costs $C_\ell$ multipole range, and the starlet its finest band

Criterion: bring the baryonic bias below 0.3$\sigma$
$C_\ell$: $\ell_{\rm max}$ tuned per survey area, from 860 at 2,000 deg² to 340 at full sky
starlet: the contamination is concentrated in the finest band, so dropping $j=1$ is enough at every area

power-spectrum baryon tension versus the upper scale cut, one panel per survey area, with the chosen cut marked against the 0.3 sigma threshold

$C_\ell$→a cut tuned to each area

starlet→whole bands only, a more conservative cut

§4 Are HOS still useful?

posteriors on baryon-safe scales at 14,000 square degrees: the power spectrum first, then peaks, then the l1-norm, all consistent with the truth

on baryon‑safe scales

starlet $\ell_1$‑norm: ×1.8 tighter than $C_\ell$ at Stage IV, ×2.6 at full sky
different degeneracy directions in the $w_0$ planes, so the HOS stay complementary to $C_\ell$
the non-Gaussian signal survives on quasi‑linear scales

and this is conservative

the whole‑band cut is not optimised, and better baryon modelling would allow smaller scales, where the gain is larger

Back to the four questions

question 1 — making the mass map Part 1

Several algorithms turn the shapes into a map of the mass distribution in the Universe. Does the choice matter for the cosmological results?

question 2 — making the mass map Part 2

Can one algorithm be accurate and fast, with error bars we can trust, for a survey the size of Euclid?

question 3 — the summary statistics Part 3

A map is too big to compare with theory, so we reduce it to a few numbers, the summary statistics. How much cosmological information do they keep, and does it take a neural network?

question 4 — the simulations Part 4

The simulations leave out astrophysics we cannot model, on the small scales. Once those scales are removed, do these statistics still constrain the parameters better than the standard analysis?

Two about the map, two about the summary statistics.

Conclusions

Q1 Does the choice of mass-mapping algorithm matter for the cosmological results? ✓Yes. The figure of merit moves by 157 %, the reconstruction error by 4 %.
Q2 Can one algorithm be accurate and fast, with error bars we can trust? ✓Yes. PnPMass is within 1 % of the fine‑tuned networks, with the smallest calibrated error bars, trained once.
Q3 How much information do the summary statistics keep, and does it take a neural network? ✓As much as an "optimal" CNN, and without training. The joint ℓ1‑norm matches the VMIM‑trained compressor.
Q4 Once the scales the simulations get wrong are removed, do the statistics still beat the standard analysis? ✓Yes. On baryon‑safe scales the ℓ1‑norm is still ×1.8 tighter than $C_\ell$ at Stage IV, ×2.6 at full sky.

Conclusions

Q1 Does the choice of mass-mapping algorithm matter for the cosmological results? ✓Yes. The figure of merit moves by 157 %, the reconstruction error by 4 %.
Q2 Can one algorithm be accurate and fast, with error bars we can trust? ✓Yes. PnPMass is within 1 % of the fine‑tuned networks, with the smallest calibrated error bars, trained once.
Q3 Can an interpretable statistic keep as much information as a neural network? ✓Yes. The joint ℓ1‑norm matches a CNN trained to be information‑optimal, and it is analytic, without any to training.
Q4 Once the scales the simulations get wrong are removed, do the statistics still beat the standard analysis? ✓Yes. On baryon‑safe scales the ℓ1‑norm is still ×1.8 tighter than $C_\ell$ at Stage IV, ×2.6 at full sky.

Conclusions

Q1 Does the choice of mass-mapping algorithm matter for the cosmological results? ✓Yes. The figure of merit moves by 157 %, the reconstruction error by 4 %.
Q2 Can one algorithm be accurate and fast, with error bars we can trust? ✓Yes. PnPMass is within 1 % of the fine‑tuned networks, with the smallest calibrated error bars, trained once.
Q3 How much information do the summary statistics keep, and does it take a neural network? ✓As much as an "optimal" CNN, and without training. The joint ℓ1‑norm matches the VMIM‑trained compressor.
Q4 Once the scales the simulations get wrong are removed, do the statistics still beat the standard analysis? ✓Yes. On baryon‑safe scales the ℓ1‑norm is still ×1.8 tighter than $C_\ell$ at Stage IV, ×2.6 at full sky.

Conclusions

Q1 Do we need deep learning to extract the non-Gaussian information? ✓No. Read the bins jointly and a fixed wavelet statistic matches the optimal compressor, with no training.
Q2 Does baryonic feedback put it out of reach? ✓No. Cut every contaminated scale and the ℓ1-norm is still ×1.8 tighter at Stage IV, ×2.6 at full sky.
Q3 Can redshift nulling then be used with higher-order statistics? ✓Yes, provided the summary reads the bins jointly. There is no loss of information, and the inflation is a frame artifact. What is left: the genuinely three- and four-bin structure a pairwise statistic cannot reach.

Backup

Cosmology

The late-time S₈ and the local H₀ disagree with the CMB predictions

compilation of S8 measurements from 2019 to 2026 against the combined CMB band — S₈ from late-time probes lies below the CMB valuePantos & Perivolaropoulos 2026, arXiv:2602.12238

compilation of Hubble constant determinations: model-dependent early-Universe values against direct distance-ladder measurements — H₀ from the early Universe against the local distance ladderDi Valentino, Levi Said & Saridakis 2025, arXiv:2509.25288

possible explanations new physics, a statistical fluctuation, or a systematic in the analysis

It was the third one

Whether a tension reflects new physics or an unrecognised systematic is settled by the analysis as much as by the data.

KiDS re-analysed with better shape and redshift calibration, and the tension fell below 1σ. And when the statistical errors shrink by an order of magnitude, that is what is left.

ΛCDM fits almost all of it — and three things in this picture are still open

history of the Universe, annotated with the open questions

Λ is a name for three different possibilities, and only precision separates them

$G_{\mu\nu}$

Modified gravitychange the left-hand side

$=$

$\dfrac{8\pi G}{c^{4}}\,T_{\mu\nu}$

Dark energya fluid with $P = w\rho$, $w \approx -1$

$-$

$\Lambda\, g_{\mu\nu}$

Cosmological constanta constant of nature

Lensing surveys measured a lower clustering amplitude than the CMB predicted

the S8 tension: DES Y3 and KiDS-1000 cosmic shear against Planck in the Omega_m S_8 plane — DES Y3 + KiDS-1000 joint reanalysis, Abbott et al. 2023

$S_8 \equiv \sigma_8\sqrt{\Omega_{\mathrm{m}}/0.3}$ — the combination cosmic shear measures best. The joint reanalysis sat 1.7σ below Planck.

Three ways to read a tension like this:

New physics
A statistical fluctuation
Systematics in the analysis

KiDS-Legacy re-analysed in 2025 and found $S_8 = 0.815$, below 1σ. It was the third one.

What do we want to know about the Universe?

The Beginning

$A_s,\, n_s,\, f_{NL},\, \tau \ldots$

Statistics of initial conditions

The Content

$\Omega_m,\, \Omega_b,\, m_\nu \ldots$

Cosmic ingredients

The Future

$H_0,\, \Omega_\Lambda,\, w_0,\, w_a \ldots$

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?
Most of it is dark, and invisible to us.

energy and matter budget of the Universe

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:

four cosmological probes: the CMB, baryon acoustic oscillations, galaxy clustering, and gravitational lensing — Clockwise: **CMB** — initial conditions; **BAO & SNe** — geometry; **Galaxy clustering** — matter growth; **Weak lensing** — total matter

The cosmological parameters, and which ones lensing sees

$\kappa$ is a projection of the matter field, so what lensing measures is how much matter there is and how clumped it is.

parameter	what it is	lensing
$\Omega_m$	the fraction of the Universe that is matter	strong
$\sigma_8$	how clumpy that matter is, on $8\,h^{-1}$Mpc scales	strong
$S_8$	$\sigma_8\sqrt{\Omega_m/0.3}$, the combination those two collapse into	what we measure
$w_0$	the dark energy equation of state, $p = w\rho$	moderate
$h$	the expansion rate today	weak
$n_s$	the tilt of the primordial power spectrum	weak
$\Omega_b$	how much of that matter is baryons	weak

The amplitude of $\kappa$ goes roughly as $\sigma_8\Omega_m^{0.5}$, so that combination is pinned and the two separately are not. Shape parameters are the CMB’s job.

only if asked — and a good place to send “why do you always quote S8?”
▲ the reason, in one line — lensing sees the amplitude of the projected matter field, and that amplitude is set by how much matter there is times how clumped it is. One number, not two
which is why the contour is a banana — along σ8 Ωm^0.5 it is tight, across it is not. S8 is just that direction given a name
dark energy — enters twice, through the distances and through how fast structure grows. Real sensitivity, weaker than the other two. Our figure of merit is in Om, sigma8, w0
if pushed on h, ns, Omega_b — these change the shape of the spectrum rather than its amplitude, and the CMB measures shape far better. We take priors on them
if asked about neutrinos — massive neutrinos suppress small-scale power, so lensing does carry some sensitivity, and higher-order statistics more than the power spectrum. Not a parameter we vary here

Born, Limber, ray tracing: three different things

Two are about the path light takes through a simulation. The third is about the projection integral in an analytic prediction.

	what it assumes	where it stands
Born	the potential is read on the straight, undeflected ray	first order in $\Phi$; post-Born terms on $\kappa$ sit well below current errors
ray tracing	nothing: rays are deflected plane by plane, along their true path	keeps the lens–lens coupling; kept for the shear and the post-Born terms
Limber	only transverse modes reach the projection, $k = \ell / f_K(\chi)$	better than a percent on Stage IV scales

Our maps are Born: CosmoGridV1 sums lens planes along the straight ray. Limber enters only where an analytic $P_\kappa$ is needed, as in the Wiener prior; every statistic here is measured on maps.

only if asked — and the confusion is worth heading off: these are not three levels of the same approximation
▲ the split — Born and ray tracing are two ways to turn a simulation into a map. Limber is an approximation in the theory prediction. Different axes
Born, in one line — evaluate the potential along the straight line instead of the bent one. First order in the potential
ray tracing, in one line — push the ray through the planes one at a time and let it bend, so it samples the mass it actually passed
why we can live with Born — the post-Born corrections to the convergence are far below the statistical errors of current surveys. They matter more for the shear, and more for higher-order statistics than for the power spectrum
if pushed on Limber — it fails at low multipoles and for narrow kernels. Neither bites here: our kernels are broad and we are not at low l

backup — not part of the talk

Mass mapping methods

Kaiser–Squires step by step, inpainting, mass mapping as an inverse problem, the Bayesian view of the prior, and MCALens.

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.

Inpainting fills the mask before the inversion, with a sparsity prior in the DCT

Kaiser–Squires inverts in Fourier space, where a gap left at zero is not a local problem: it spreads over every mode. Border artefacts and E/B leakage follow.

The idea

A mask costs you sparsity: its cut edges are badly represented in the DCT and throw off spurious coefficients. So take the map with the fewest coefficients that still matches the data outside the mask.

$$\min_{\kappa}\; \lVert \Phi^{\mathsf{T}}\kappa \rVert_0 \quad \text{s.t.} \quad \lVert \kappa_{\rm obs} - \mathbf{M}\kappa \rVert_2 \le \sigma$$

Solved by iterative thresholding

$$\kappa^{\,n+1} \;=\; \Delta_{\Phi,\lambda_n}\!\left( \kappa^{\,n} + \mathbf{M}\big(\kappa_{\rm obs} - \kappa^{\,n}\big)\right)$$

Decompose, hard-threshold at $\lambda_n$, reconstruct, put the observed pixels back. Repeat. $\lambda_n$ starts at the largest coefficient and cools to zero over about a hundred iterations, so the strong coefficients go in first and the faint detail last.

the same field reconstructed by Kaiser-Squires and by iterative Kaiser-Squires; in the KS panel the masked regions are flat, and after inpainting they carry structure matching the rest of the field

Why the DCT

Pires et al. tried wavelets, curvelets and their combinations. Wavelets fit an isolated halo better. But most of a κ map is the texture between the haloes, and there the DCT wins.

Pires et al. 2009 · Elad et al. 2005

Inpainting is preprocessing, so it bolts onto any method: inpainted Wiener, sparse recovery and MCALens variants all exist. A method with its own prior can instead put $\mathbf{M}$ in the forward operator and let the reconstruction fill the gap.

why it is needed — KS goes through Fourier, so a hole left at zero contaminates every scale at once, not just its own patch.
the idea, in words — a map with a hole cut in it needs far more DCT coefficients than the same map without one. Ask for the filling that needs the fewest, and you get one that looks like the rest of the field.
the iteration, in words — transform, keep the coefficients above the threshold, transform back, restore the real pixels, lower the threshold, repeat. Big structures land first, fine detail last.
the DCT — they tried wavelets and curvelets too. Wavelets are better on a single halo, but a κ map is mostly the texture in between.

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.

Every mass-mapping method is a Bayesian inference with a different prior on κ

The posterior probability of a map given the shear is proportional to the likelihood of the shear given the map, times the prior on the map.

after the data

$p(\kappa \mid \gamma)$

posterior

∝

how well the map explains the shear

$p(\gamma \mid \kappa)$

likelihood (forward model)

before the data

$p(\kappa)$

prior

1Wiener filter κ is a Gaussian random field

2sparse recovery κ is sparse in a wavelet basis

3MCALens κ is Gaussian plus non‑Gaussian

4deep learning prior learned from simulations

The same structure returns in Part 3, with the cosmological parameters as the unknown and no explicit likelihood.

The Bayesian view makes the prior explicit

$p(\kappa \mid \gamma)$

posterior

the map, given the data

what we want

∝

$p(\gamma \mid \kappa)$

likelihood

how well a map explains the measured shear

the forward model

$p(\kappa)$

prior

what we assume about the map

the only term that differs between methods

Wiener filter

κ is a Gaussian random field

sparse recovery

κ is sparse in a wavelet basis

MCALens

κ is Gaussian plus non‑Gaussian

deep learning

the prior itself is learned from simulations

The same structure returns in Part 3, with the cosmological parameters as the unknown and no explicit likelihood.

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

Bayesian reconstruction

Mass mapping problem → statistical inference problem
Goal: infer most probable value of $\kappa$-field given observed shear data

$$ \underbrace{p(\boldsymbol{\kappa} \mid \boldsymbol{\gamma}, \mathcal{M})}_{\text {Posterior }} \propto \underbrace{p(\gamma \mid \boldsymbol{\kappa}, \mathcal{M})}_{\text {likelihood }} \,\, \underbrace{p(\boldsymbol{\kappa} \mid \mathcal{M})}_{\text {prior }} $$ $\mathcal{M}$: cosmological model

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.

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 — the algebra

Models $\kappa$-field as a sum of a Gaussian and non-Gaussian component $$\boxed{\kappa = \underbrace{\kappa_{\rm G}}_{\text{Gaussian: Wiener filter}} + \underbrace{\kappa_{\rm NG}}_{\text{non-Gaussian: sparse wavelet}}}$$ $$\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)$$

backup — not part of the talk

PnPMass

The denoiser and the transformer it is built from, the residual variant, the per-pixel uncertainties and their conformal calibration, and the implementation.

The denoiser is trained once, on white noise, and never sees the mask or the covariance

The forward step carries the physics. The denoiser only has to know what a plausible convergence field looks like.

The network

SUNet: a U-Net whose blocks are Swin-Transformer blocks rather than convolutions, with 7.2 M parameters. The noise level σ goes in as an extra input channel, so one model covers a range of noise instead of one level.

a U-Net: an encoder contracting to a bottleneck, a decoder expanding back to full resolution, and skip connections carrying fine detail across at each scale

How it learns

Take a simulated κ, add white Gaussian noise, ask for the κ back.

$$\hat\theta \;=\; \arg\min_{\theta}\; \frac{1}{n} \sum_i \big\lVert D_{\theta}\big(\kappa_i + n_i\big) - \kappa_i \big\rVert^2$$

σ is drawn uniformly over 0 to 0.2. No shear, no mask and no survey enter the training set.

Why the mask and the covariance drop out

The forward step may use any linear operator. PnPMass takes $\mathbf{B} = \mathbf{A}^{\!\top}\Sigma^{-1/2}$, which whitens the residual before adding it back. The denoiser then sees κ plus white noise of width $\tau$, whatever mask and galaxy density went into $\Sigma$. A MAP step would use $\Sigma^{-1}$, and the noise would arrive coloured.

Leterme, Tersenov, Fadili & Starck 2026, A&A 710, A292

What survives is the scalar $\tau$, whose usable range depends on $\Sigma$. DeepMass needs the full covariance, and a retraining for every footprint.

its job — the loop alternates a data step and this. The data step knows about shear, mask and noise; the denoiser only knows what a mass map looks like.
how it learns, in words — nothing clever. Show it a simulated map with noise added, ask for the map back, repeat a few hundred thousand times. It is not trained on the lensing problem at all.
the trick — you get to choose the operator in the forward step. Choosing the minus one half power of the covariance whitens the residual, so the noise arriving at the denoiser is white whatever the survey did. Half power, not the inverse the likelihood would give you.
if pushed — one scalar is left, the step size, and the network takes the noise level as an input, so a single training covers the range.

A transformer lets every patch look at every other patch

top: the local receptive field of a convolution against the global receptive field of self-attention, drawn on a convergence map; bottom: one self-attention update, where a query patch is compared with all patches to give attention weights, and the output is their weighted sum

One update

Cut the map into patches. Compare one patch with all the others, turn the comparisons into weights, and replace it with the weighted average of the rest.

Why it is worth the trouble

The weights come from the content of the patches, not from where they sit, so the network learns which distant regions belong together. A convolution can only ever look at a fixed neighbourhood.

The catch, and the fix

Comparing every pair costs time growing as the square of the number of patches. The Swin variant keeps attention inside local windows that shift between layers, so global reach is rebuilt over several layers.

Vaswani et al. 2017 · Liu et al. 2021

PnPMass on residuals

§2 Pixel-wise uncertainties from a second network, calibrated with conformal prediction

shear $\gamma$the data

$\hat\kappa = F_{\Theta}(\gamma)$the PnPMass map

$\hat\sigma^{2} = G_{\Omega}(\gamma)$a second networkthis work

calibrated intervalCQR, held‑out set

where the spread comes from The posterior spread comes from the noise and the mask (Part 1), and from the learned prior that selects one map. The predicted variance has to account for both.

first — predict the error

Train $G_{\Omega}$ on the squared residual of the reconstruction,

$$\min_{\Omega}\ \mathbb{E}\left[\left\| G_{\Omega}(\gamma) - \big(\kappa - F_{\Theta}(\gamma)\big)^{2} \right\|_2^{2}\right]$$

whose minimiser is the posterior variance, pixel by pixel. One forward pass, no sampling.

then — calibrate it

The variance predicted by a network is not calibrated.

Conformal quantile regression rescales it pixel by pixel on a held‑out calibration set, so that the intervals reach the target coverage.

The result is a coverage guarantee that holds whether or not the network is well specified.

the chain — shear in, the denoiser gives the map, a second network trained on the same simulated pairs gives the error
» where the spread comes from — noise and mask leave a whole family of maps consistent with the data — the fan from Part 1 — and the learned prior picks one of them. Both have to be covered
how — trained on the squared residual of the reconstruction. The minimiser of that loss is the posterior variance. One forward pass, no sampling
» then calibrate — a network's own variance is overconfident. CQR, a held-out set, a per-pixel rescaling
» ▲ the shift — not “the network says plus or minus sigma”, but a stated coverage level that holds anyway

Conformal calibration fixes the size of the error bar, without trusting the network

A network's error bar comes out the wrong size, and nothing tells you by how much. This is how you fix it.

three panels: the network's interval drawn on fourteen held-out examples with seven truths landing outside; the sorted conformity scores with the quantile q marked; the same intervals widened by q, with three truths outside, the requested rate

What you get

$$\alpha - \tfrac{1}{n+1} \;\le\; \mathbb{P}\big\{\kappa \notin [\hat\kappa^-,\, \hat\kappa^+]\big\} \;\le\; \alpha$$

That fraction falls outside, whatever produced the bar. The calibration set can be small, nothing has to be Gaussian, and the network does not have to be right. The one condition is exchangeability: calibration maps and real data drawn on equal terms.

What you do not get

The guarantee is marginal: an average over many maps, not a promise about the one in front of you. Some regions do worse than the average, and peaks are the known case — the Wiener bars miss there systematically even after calibration.

Here: $2\sigma$, so $\alpha \approx 4.55\%$, from 1024 calibration maps. The correction sets a minimum size for that set: 21 maps at $2\sigma$, 370 at $3\sigma$, 15 787 at $4\sigma$.

Romano et al. 2019 · Leterme, Tersenov, Fadili & Starck 2026

the problem, in words — the second network predicts a variance, and you could call that a 2σ bar. But nothing forces a network's variance to be right, so the bar is the wrong size and you do not know by how much.
the fix, in words — hold out maps where you know the truth. For each one, ask how far outside the bar the truth landed. Take the 95th percentile of that, and widen every bar by it. Now 5% land outside, by construction.
why it is strong — you never assumed the network was right, or that anything was Gaussian. A terrible network gives honest but very wide bars.
it can also shrink — if the bars were too conservative the quantile is negative. That is what happened to DeepMass, whose bars came out smaller after calibration.
the honest limit — marginal, not conditional. On average over maps, not guaranteed on the one map you have, and peaks are where it shows.

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

backup — not part of the talk

Method detail

Kaiser–Squires step by step, the Bayesian reconstruction stack, sparse recovery, the MCALens algebra, and the PnPMass implementation. Here for questions.

Wavelets are localised in both scale and position

left: three sines, each one frequency spanning the whole axis; right: the starlet wavelet at three dilations, compact and localised

a wavelet is a localised, oscillating function $\psi$ with zero mean, $\int\psi\,\mathrm{d}x = 0$
dilations and translations form a family, $\psi_{a,b}(x) = \tfrac{1}{\sqrt{a}}\, \psi\!\left(\tfrac{x-b}{a}\right)$; the coefficient $W(a,b) = \int f\,\psi_{a,b}$ measures the structure at scale $a$ and position $b$
a sine has infinite support: full frequency resolution, no spatial resolution. A wavelet trades some of the former for the latter
for a field of localised structures (peaks, filaments, voids) this is the right trade. The power spectrum is the Fourier statistic; peaks and the ℓ₁‑norm are the wavelet ones

only if asked
start from Fourier — sines are one exact frequency spread over the whole domain, so a coefficient says which scales are present and nothing about where
a wavelet is the other trade — a small wave: it oscillates, then dies away, and has zero mean
the family — dilate and shift one shape. The coefficient is the overlap, so it says how much structure of that size sits at that place
zero mean is what makes it a detector — over a flat region every coefficient is zero
the price — a little frequency resolution. For a field of localised objects that is a bargain
▲ the one-line summary — the power spectrum is the Fourier answer, peaks and the l1-norm are the wavelet one

The starlet transform splits the data into one band-pass image per scale plus a coarse residual

a signal made of three wave packets, decomposed into five starlet bands plus a coarse residual; each packet appears in the band matching its size

convolving the data with the dilated wavelet at each scale gives one band‑pass image per scale, plus a coarse residual
the starlet is isotropic and undecimated: every band keeps the full pixel grid, and the reconstruction is exact, $I = c_J + \sum_j w_j$
a structure of a given angular size appears in a single band, so a scale cut amounts to dropping one band

In image processing, the same transform is used for compression and denoising

three levels of a 2-D wavelet transform of a photograph, in the classical quadrant mosaic: the approximation shrinking into the corner, with vertical, horizontal and diagonal detail bands at each level

In two dimensions each level splits into vertical, horizontal and diagonal details, and recurses on what is left.

compression: most of the image is in a few large coefficients; keep those and discard the rest (JPEG 2000)
denoising: threshold the small coefficients (the proximal step of Part 1)
the starlet is isotropic and undecimated: no preferred direction, and every scale on the full pixel grid

Why filter at all, and which filter?

Core motivations

Pixel domain mixes scales → hard to attribute physics by scale.
Filtering = change of basis → separates spatial frequencies; reduces mode mixing.
Non-stationarity: WL $\kappa$ is scale-variant & non-Gaussian → need localization in space+scale.
Sparsity: the right basis concentrates energy in few coefficients → higher SNR for HOS tails.

Comparative view

Fourier: perfect spectral selectivity, but global (no spatial localization).
Gaussian: low-pass; simple, but mixes nearby scales and blurs features.
Starlet (à trous): isotropic, undecimated, shift-invariant multiscale bands.
Update rules: $c_{j+1} = c_j * h_j$, $w_{j+1} = c_j - c_{j+1}$; reconstruct with $\sum_j w_j + c_J$.

The bands are band-passes and adjacent ones overlap in multipole, so the covariance comes out nearly diagonal rather than diagonal.

The non-Gaussian information is in the phases

the original density field, the same field keeping only the phases, and keeping only the amplitudes

the power spectrum keeps only the Fourier amplitudes
the phases carry the structure of the cosmic web (filaments, haloes, voids)
higher-order statistics are needed to access it: peak counts, wavelet statistics, the ℓ₁‑norm, Minkowski functionals

Peak counts: the maxima of the signal‑to‑noise field

a convergence map read by the power spectrum and by peak counts

A peak is a local maximum of $\nu = \dfrac{\left(\mathcal{W} \ast \kappa\right)(\theta_{\rm ker})}{\sigma_n^{\rm filt}}$ , the signal-to-noise field of the filtered map.

peaks trace regions of high $\kappa$, i.e. massive structures
the peak function (counts per SNR bin) is a one‑point statistic sensitive to non-Gaussian information

A single starlet transform gives peak counts at every scale

a map written as a sum of starlet band-pass images plus a coarse map — The starlet transform writes a map as a sum of band‑pass images, each carrying structure of one characteristic angular size, plus a coarse map.

counting peaks band by band gives a multi‑scale analysis
all scales from a single transform
the bands cover different frequency ranges, so the peak‑count covariance is almost diagonal

The starlet $\ell_1$‑norm uses every pixel

the starlet decomposition of a weak lensing map, band by band

$$\ell_1^{j,i} = \sum_{u=1}^{\textrm{\# coef}(\mathcal{S}_{i,j})} \left| \mathcal{S}_{j,i}[u] \right| = \| \mathcal{S}_{j,i} \|_1$$

the sum of absolute starlet coefficients, per band and per SNR bin
peaks and voids are included automatically, with no threshold and no peak definition
information from all the pixels, at every scale

Asymmetry of $\delta$ and its imprint on wavelet $\ell_1$

$$\delta = \frac{\rho - \bar\rho}{\bar\rho}, \qquad \delta \ge -1 \;\; (\rho \to 0), \qquad \delta \to +\infty \;\; (\rho \gg \bar\rho)$$

Gravitational collapse produces unbounded high-density peaks, but voids cannot go below −1.
This asymmetry yields larger wavelet amplitudes in the positive tail than in the negative tail.
Thus the wavelet $\ell_1$-norm at a given scale is dominated by rare collapsed peaks; unlike peak counts it still reads the underdense side, in its negative SNR bins.
$\kappa$ is a projection of $\delta$ → bounded below at the empty beam, unbounded above, and Gaussianized by the sum along the line of sight.

Why $\ell_1$ instead of $\ell_2$ on wavelet coefficients?

Robust to tails: $\ell_2$ squares values → the sum is driven by rare extremes; $\ell_1$ grows linearly → the tail contributes in proportion.
Sparsity: non-Gaussian structure is sparse in wavelets; $\ell_1$ concentrates information, $\ell_2$ dilutes it across many small modes.
Binning does the work: $\ell_1$ is summed per band and per SNR bin → the bins trace the shape of the coefficient distribution.
Cost vs alignment: $\ell_1$ is non-smooth and harder to optimise, but its geometry matches the sparsity of the field we measure.

→ $\ell_1$ maximizes cosmological signal-to-noise on non-Gaussian structure.

Ajani, Starck & Pettorino 2021, the paper that defines the starlet ℓ₁‑norm.

PDF vs wavelet $\ell_1$: differentiability matters

Key distinction

Histogram / PDF (measured) is a piecewise-constant functional of samples → zero or undefined gradients at bin edges.
Wavelet $\ell_1 = \sum_i |w_i|$ is sub-differentiable everywhere.
Thus $\ell_1$ can flow gradients in any differentiable forward / inverse pipeline; PDFs cannot without a kernel or a soft-binning reparametrization.
Interpretation: the PDF is a descriptive statistic; $\ell_1$ is also an optimizable functional of the field.

Bin edges are still a choice in both cases, and our derivatives with respect to cosmology are finite differences on simulations either way.

Does filtering bias inference?

Inference integrity

Filtering is part of the summary definition, not a treatment applied to the data → the same transform runs on the observed map and on every mock.
It biases only when applied on one side: to the data but not the prediction, or with a mask and apodisation that differ between data and mocks.
Covariance: estimated from a finite set of realizations → a noisy $\mathbf{C}^{-1}$ moves the contours, whatever the filter.
Likelihood: no analytic mean for peaks or $\ell_1$ → simulation-based inference rather than an assumed Gaussian.
Diagnostics: TARP-DRP coverage and SBC rank uniformity, every arm, in the standard frame and the nulled one.

Which makes it a measurement rather than an argument. The coverage card in Part 2 is the answer.

Physical scales in cosmology

	comoving size [$h^{-1}$Mpc]
angular scale	$z=0.1$	$z=0.3$	$z=0.6$
1′	0.09	0.25	0.46
10′	0.9	2.4	4.6
20′	1.7	4.9	9.1
1°	5.1	15	27
structure
galaxy halo	0.1 – 0.3
cluster, $R_{200}$	1 – 2
$\sigma_8$ sphere, non-linear scale	8
linear regime	≳ 63
BAO, sound horizon $r_d$	99

Comoving, for $\Omega_m = 0.26$, $h = 0.674$. The BAO sound horizon is 147 Mpc, which is the 99 above once divided by $h$.

the reference card — for “what is that in Mpc?”. Find the row, do not read the table
▲ say the caveat first — an angle is a size only once you fix the redshift. That is why there are three columns
the one-line answer — arcminute scales are sub-Mpc to a few Mpc. Our bands, 10 to 80 arcmin, run from about 1 to 40
the anchor — 8 h⁻¹ Mpc is the sigma-8 sphere and the turn to non-linear. That is the scale everyone in the room already has
if asked about baryons — feedback bites below about 1 h⁻¹ Mpc, which is the finest band and nothing else. That is why one band is the whole scale cut
if asked about BAO — 147 Mpc is the number people remember, and it is 99 in h⁻¹ Mpc. Either way far above anything we measure

Part 2 of 2

Learned vs analytical,
and can we trust it?

The analytical ℓ1-norm vs a learned CNN, calibration, and the answer to the BNT puzzle.

Q1 The joint ℓ1-norm reaches the FoM of the optimal compressor

three-parameter figure of merit, matched pipeline

the three-parameter figure of merit for four summaries, drawn as stems with error bars: l1 on auto-maps 2448, l1 plus product cross-maps 3045, the joint l1-norm 3371 and the VMIM-trained CNN 3326, the last two agreeing within their errors

per-mock distributions of the marginal uncertainties and the figure of merit for the l1 plus product, joint l1 and CNN summaries; the medians coincide — a tie on every parameter over 9,000 mock observations, so not an artefact of the figure of merit

Q1 At equal constraining power, the analytical statistic is cheaper and safer to use

constraining power joint ℓ1-norm 3371 ≈ CNN 3326

analytical — the ℓ1-norm

Fixed before any simulation is seen
No training, nothing to overtrain
Survives a change of survey or forward model
Open to inspection, scale by scale
Analytically predictable on large scales

learned — the compressor

Retrained for every change of configuration
Seed-to-seed scatter
Ten coordinates with no meaning to inspect
Cannot tell physics from simulator artefacts

One advantage of the compressor remains: it reads the bins jointly by construction. This matters in the nulled frame.

Learning optimal filters vs fixed multiscale starlet

Learned summaries in lensing:
– Gupta et al. (2018), CNN on noiseless $\kappa$ maps, $\times 5$ over $C_\ell$ and $\times 4$ over peaks PRD 97, 103515
– Fluri et al. (2018), the same with shape noise PRD 98, 123518
– Ribli et al. (2019), beats peak counts, then points back at an analytic statistic: peak steepness Nat. Astron. 3, 93
– Jeffrey, Alsing & Lanusse (2021), neural compression for likelihood-free inference on DES SV MNRAS 501, 954
Trade-offs of learned filters:
– Address one task but become training-dataset dependent.
– Lose scale-by-scale interpretability and analytic covariance modelling.
– No formula to publish; reproduction needs the weights.
Why we retain the starlet:
– Isotropic, undecimated multiresolution bands (clear physics per scale).
– Fast algorithmic complexity $\mathcal{O}(N J)$, band responses known.
– Amenable to analytic summary modelling and systematic control for pipelines like Euclid.

backup We are all optimizing statistics; the two-point camp still does not trust the contours

🧐

the 2-point camp

"I don't believe any of your contours."

what would make HOS flagship-grade?

blinding
robust covariance
emulators
systematics
analytical cross-checks
non-Gaussian likelihood
method limits
null / validation tests
simplicity

§3 Part 2: learned summaries, and the BNT cliffhanger

the question How much better are "optimal", learned summaries than our hand-built summary statistics?

what is a learned summary?

A neural network that compresses the κ map directly into a few numbers, instead of a hand-designed statistic
Trained with VMIM to keep the cosmological information: the "optimal learned compressor"

and, left over from Part 1 ...and what the hell is going on with BNT?

Training a neural summary I: regression (MSE)

simulator

$\theta\sim$ prior

→

map $x$

→

network

$f_\phi$

→

estimate

$\hat\theta$

\[\mathcal{L}=\mathbb{E}\,\big\lVert\,\theta-f_\phi(x)\,\big\rVert^{2}\]

parameter space $\theta$

§3 Can we trust it? Every arm passes the same TARP + SBC tests reasonably

TARP-DRP expected coverage against credibility level, every arm arriving in turn and every one lying on the diagonal — TARP-DRP coverage: on the diagonal

SBC rank densities for every arm and every parameter, flat within the band, in both the standard and the nulled frame — SBC ranks: flat within the band, standard frame and nulled (dashed)

SBC ranks: flat within the band, standard frame and nulled (dashed)

tight is not the same as correct

Both arms pass the same battery: varied-θ TARP-DRP coverage + SBC rank uniformity
The constraining power is real

TARP and SBC: do the error bars mean what they say?

When the pipeline says it is 68% confident, is it right 68% of the time? Both tests answer that by rerunning the analysis on simulations where the true parameters are known.

SBC, one parameter at a time

For each simulation, count how many of the pipeline's sampled values fall below the true one.
If the pipeline is honest the true value is just one more draw from what it produced, so that count is equally likely to be anything.
Counts piling up at the extremes mean the error bars are too small; piling up in the middle means too large.
Talts et al. 2018

TARP, all of them at once

Runs the same check on all parameters together, using distances to a randomly chosen point in parameter space.
Comes with a proof: a posterior that passes TARP is correct. Passing the one-parameter test does not guarantee that.
Lemos et al. 2023

Both average over many simulated datasets: they vouch for the pipeline in general, not for the one dataset we have.

only if asked — or if someone challenges the contours rather than the method
▲ the question both answer — when we say 68%, are we right 68% of the time
SBC, the mechanism — that count is the rank. A correct posterior makes the truth just another draw from it, so the rank histogram over many simulations comes out flat. U-shaped means too narrow, humped means too wide, sloped means biased
TARP, the mechanism — pick a random point, count the samples closer to it than the truth is. That fraction is the smallest credible level that would have caught the truth, and averaging it gives the coverage actually achieved
why both — they fail differently. SBC on every parameter separately can pass while the joint posterior is still wrong, which is what the TARP proof closes
the caveat, say it before they do — both average over simulations, so they certify the pipeline, not this dataset. L-C2ST is the local version where the dimension allows it

§3 The analytical statistic ties the CNN, at lower cost and risk

the benefit

a tie

Joint ℓ1 3371, CNN 3326, both calibrated. The optimised network buys no measurable constraining power.

the cost

extensive architecture and hyperparameter search
a very large dataset (899 cosmologies, 3.2×10⁵ patches); VMIM needs the scale or it biases
unphysical-information traps (patch geometry, map-mean / mass-sheet mode, 20° projection features) that tighten contours dishonestly
and these largely escape TARP and SBC: the contours look calibrated and are still wrong

in favour of the ℓ1-norm

ℓ1 is simple, interpretable, inspectable; CNNs are powerful but treacherous
Where the CNN earns its keep: BNT, the channel-mixing win
So the question is what the cost and the risk are buying

§3 Cross-bin information: the κ_iκ_j product gives +24%, the joint ℓ1-norm the rest

the completeness ladder in figure of merit: auto 2448, plus convolution 2671, plus product 3045, both 3255, joint l1 3371

cross-map constructions

Product κ_iκ_j (= ξ_ij): +24% over auto-only
A convolution buys +9%, and is sensitive to the training realisation
2448 → 2671 → 3045 → 3255, and the joint ℓ1 closes it at 3371

A full-sphere cross construction inflates this, but its channels carry 12–20% of their variance at ℓ < 18 against 0.4–1% for the autos. That is leakage. We use only the physically buildable flat-sky arms.

backup Where the ℓ1-norm’s constraining power comes from

relative Fisher information of the auto starlet l1-norm per wavelet scale and signal-to-noise bin, for each of the three parameters

Intermediate scales ($j = 2$–3, tens of arcminutes) and moderate signal-to-noise ($|\nu| \approx 1$–2): mildly non-linear structure rather than the rare extreme peaks. Consistent with Paper I: drop $j=1$ and you keep most of it.

backup The ℓ1-norm across cosmologies

starlet l1-norm data vectors coloured by sigma-8, per wavelet scale and tomographic bin

Coloured by $\sigma_8$. The statistic responds smoothly and monotonically across the prior: no training, no emulator, and the dependence is visible by eye.

backup The data: flat-sky tomographic patches

example 10 degree gnomonic convergence patches across the four tomographic bins, noiseless and noisy

10°×10° gnomonic patches, 80×80 px (7.5′/px), $|b| < 75^{\circ}$. 180 patches per realisation × 50 noise realisations = 9,000 mock observations. Patch size chosen at 10° because gnomonic corner distortion falls from 6.3% at 20° to 1.5% at 10°.

backup The tie per mock observation

per-mock distributions of the three marginal uncertainties and the figure of merit, for the l1 plus product, joint l1 and CNN arms

σ(Ω_m), σ(σ₈), σ(w₀) and FoM₃ across 9,000 mocks: medians 3045, 3371, 3326

The answer to “FoM₃ is fragile in a degenerate parameter space”: the marginals are quoted alongside it, and this is the population rather than a point. The ± on the medians is the spread over three independently trained compressors, the dominant source of run-to-run variability.

Part 1 of 2

Do baryons break HOS?

Baryonic feedback, the wavelet ℓ1-norm, and the BNT transform.

Q2 Without any feedback model, the ℓ1-norm stays ahead of the power spectrum

figure of merit relative to the power spectrum, on each statistic's own baryon-safe scales

Conservative case: no feedback model, every contaminated scale removed.

§4 Stage IV is no longer statistics-limited, it is systematics-limited

to trust a statistic

Before we trust any summary statistic, we have to quantify how each systematic affects it, and at the contour level (the inferred parameters).

Illustris: baryonic feedback reshaping the cosmic web

the systematic at hand

Baryonic feedback (AGN, supernovae) suppresses matter on small scales, mimicking cosmological signal and biasing inference, exactly where the constraining power lives and where the feedback models disagree most.

core questions

How does unmodeled baryonic feedback bias our non-Gaussian statistics?
After safe scale cuts, do HOS still outperform the power spectrum?

§4 Nulling promises localized scale cuts... For HOS it inflates the contours instead

the promise

standard lensing efficiency kernels: broad and overlapping across redshift — broad, overlapping

↓

BNT lensing efficiency kernels: nulled and localized in redshift — nulled, localized

A fixed angular scale stops mixing redshifts, so the cut can go only where the systematic is

the problem

the l1-norm posterior in the BNT basis, visibly inflated against the same statistic in the standard basis — blue = the ℓ1-norm in the nulled basis

BNT = a fixed, invertible transform → nothing can have been lost...

Tersenov+ 2026, A&A and Tersenov+, submitted

what BNT is — Bernardeau, Nishimichi & Taruya 2014. A linear re-mixing of the tomographic bins that cancels the low-redshift lensing efficiency
why anyone wants it — standard kernels are broad and overlapping, so one angular scale mixes low-z small scales with high-z large ones, and an angular cut throws away clean high-z information along with the contamination
nulled — each transformed field is localised in redshift, so the baryon sensitivity lands in one bin and you cut scales only there
» it backfires — the same mixing correlates the originally independent shape noise across bins. The noise floor rises and the contours inflate, worse than standard tomography even with conservative cuts
plant, do not explain yet — for the power spectrum it is lossless, provided you keep the cross-spectra. It pays off on the next slide
then the paradox, and pause — this is a fixed invertible matrix, so nothing can have been lost

§4 The information is recovered by a summary that reads the bins jointly

figure of merit retained under the nulling

one bin at a timeℓ1, auto-maps

0.16×

+ one derived field per pairℓ1 + product cross-maps

0.24×

each pair's full 2-D distributionjoint ℓ1-norm

0.72×

all four channels at onceCNN, VMIM-trained

~1×

Same four summaries. In the standard frame they span 38% in FoM; in the nulled frame, a factor of six.

nulled frame over standard

posteriors in the nulled frame: the auto-only l1 balloons, retaining 0.16 of its figure of merit, then the product and joint arms tighten — dashed = the nulled frame

Nulling costs no constraining power provided some stage of the pipeline reads the bins jointly. The power spectrum with its cross-spectra is the simplest example.

the ladder — per-bin l1 keeps 0.16. Add the product cross-maps, 0.24. The joint l1, built two slides ago, keeps 0.72
the CNN — reads all four channels natively and is essentially untouched: 0.96, with an error consistent with 1
▲ how to quote it — “no measurable loss”. Give the number only if asked
▲ the practical message — nulling need not cost a higher-order analysis anything, provided some stage of the pipeline reads the bins jointly
cash the clause from the last slide — the power spectrum survives for exactly this reason. Auto plus cross spectra are closed under the transform; the autos alone are not
worth saying out loud — nulling is designed to de-correlate the bins, so joint reading matters more there, not less
the residual — 0.72 is not all of it. The remainder is genuine three- and four-bin structure a pairwise statistic cannot reach, and the one place the network keeps an advantage

Same information, a different frame: why BNT collapses the ℓ1-norm but not the CNN

the point cloud (all the information) never moves

lensing kernels q(z)

Signal & noise under BNT: who can read it

noise covariance

What survives BNT: the 2-point rule

backup Do baryons break HOS? No.

Part 1, baryons Usable non-Gaussian information persists on baryon-safe scales (the ℓ1-norm beats P(k) ×1.8 at Stage IV, ×2.6 at full sky), cleaned by a single scale cut.

Part 2, learned vs analytical Read the bins jointly and the hand-built ℓ1-norm matches the optimal learned summary: 3371 against 3326, a tie, both calibrated.

BNT The apparent BNT break is a frame artifact: what survives tracks how jointly the summary reads the bins: 0.16, 0.24, 0.72, 0.96.

the bigger question And can we trust higher-order weak lensing? Getting there...

§4 BNT amplifies and correlates the shape noise in the maps

noisy tomographic maps before BNT — before BNT

noisy maps after BNT, visibly degraded SNR — after BNT: the SNR collapses

§4 Without noise, BNT redistributes the signal without loss

noiseless tomographic maps before BNT — before BNT (noiseless)

noiseless maps after BNT — after BNT (noiseless)

Without shape noise, BNT is an invertible redistribution of the signal (the common mode becomes one shallow map plus thin slices). The contour inflation comes from the correlated noise, not from lost signal.

backup Summary: BNT becomes viable once the summary statistic mixes the bins

per-bin l1 contours inflate under BNT — problem: the per-bin ℓ1 inflates under BNT

→

the channel-mixing CNN posterior is indistinguishable in the standard and the nulled frame — resolution: a summary that reads the bins jointly is unaffected

the common thread: cross-bin information

P(k) → ℓ1 (much more, even on safe scales) → learned (a tie, both calibrated)
Per-bin statistics cannot access cross-bin info (break under BNT); a channel-mixing compressor can (BNT-lossless)
BNT becomes viable once the summary mixes bins

Outlook: a route to baryon-robust non-Gaussian SBI that keeps BNT's per-bin scale cuts without the contour inflation. A next step rather than a finished end-to-end measurement.

backup The starlet band responses in multipole

normalised squared response of each starlet band against multipole, the four dyadic bands plus the coarse scale

This is what "drop $j=1$" removes: the highest-frequency band, everything above $\ell \approx 500$. The other bands are untouched, at every survey area.

backup BNT on the power spectrum: the bin-specific cut is worth ×1.4

power-spectrum posteriors at 14,000 square degrees: the standard basis with a global lmax of 460 against the BNT basis cut only in the first transformed bin

Baryon sensitivity localises to the first transformed bin, so the cut goes there alone and bins 2–4 keep $\ell_{\rm max} \approx 1024$: 92 of 120 bandpowers retained against 50 for the global cut. σ(Ω_m) −14%, σ(σ₈) −19%.

backup Baryonic suppression of the tomographic power spectra

$fractional change in the tomographic angular power spectra from baryonification, one curve per source bin$

Reaches ≈1.5% by $\ell \approx 1000$. The apparent "high-z bins are worse" ordering is a noise artefact: the noise power dominates the denominator at low z. On noiseless maps the ordering inverts to the physically expected one.

backup Baryonic response of the ℓ1-norm, scale by scale

$fractional change in the starlet l1-norm from baryonification, per wavelet scale and signal-to-noise bin, one curve per source bin$

Concentrated in the positive SNR tail, and vanishing in the noise-dominated bulk ($|\nu| \lesssim 2.5$). That is why an SNR-space cut is available to higher-order statistics and not to the power spectrum. Left to future work in the paper.

backup Coverage in the nulled frame: every arm still passes

TARP-DRP coverage in the BNT basis, every arm on the diagonal

The collapse is calibrated: the wide contours report a real loss in that representation, not over-confidence. So the retention ladder measures information, not the quality of a fit.

§4 The channel-mixing CNN is unaffected by the transform (0.96×)

the CNN posterior in the standard and the nulled frame, the two contours lying on top of each other

same maps, same transform, no cost Feeding the network the nulled maps gives the same first layer with kernels KB, so “undo the nulling” is one configuration of the first layer, available before any non-linearity at no capacity cost. The 4% shortfall is an optimisation residual rather than lost information.

The other systematics, and where each one stands

Baryonic feedback is Part 4. These are the rest.

systematic	what it is	how it is handled
intrinsic alignments	galaxies align with the tidal field, so shapes correlate without any lensing	2pt: NLA, TATT. HOS: open
photometric redshifts	$n(z)$ from colours, not spectra. Stage IV needs $\bar z$ per bin to 0.002	shift and width nuisances
shape measurement	PSF modelling, multiplicative and additive shear calibration, blending	calibrated on image simulations
masks and geometry	Kaiser–Squires is non-local, so gaps leak across the map: spurious peaks and B‑modes	inpainting, regularised mapping
covariance	no closed form for HOS. Hartlap factor, super-sample covariance	large simulation suites
source clustering	sources are assumed unclustered, and uncorrelated with the lenses	sub-percent for 2pt

Every one of these is better understood for the two-point function than for higher-order statistics, and all of them bite hardest on the small scales where HOS earn their advantage.

only if asked — “what about the other systematics?”
▲ the honest summary — baryons are the one we took on. The others are real, and for higher-order statistics most of them are less well characterised than for the power spectrum
intrinsic alignments, if pressed — two terms. II correlates neighbours in the same tidal field and tomography suppresses it. GI correlates a foreground shape with the shear it induces behind, survives redshift separation, and is degenerate with the amplitude of structure. That is the one that matters
why IA is the open one for us — the alignment signal is itself non-Gaussian, so a nuisance amplitude fitted to the 2-point function does not tell you how it propagates into peaks or the l1-norm
masks are our own subject — that is what the inpainting card in the mass-mapping column is about
if they ask which worries you most — intrinsic alignments, for the reason above. Say so plainly

What comes next

All of it on simulations, deliberately. That is the first thing to change.

Near term

To data — UNIONS, Euclid. Every systematic built into the forward model at map level.
Close the loop — PnPMass maps and their pixel uncertainties → a higher-order inference, end to end.
Learned maps, same pipeline — what survives away from the training cosmology. Fix by training across them, or conditioning the denoiser.

Further out

PnPMass, next versions — spherical for Stage IV, tomographic for the inter-bin correlations, BNT-aware for cross-bin information in the map.
Conditional coverage — the conformal guarantee where the miscoverage sits, at the peaks.
The systematics left out — alignments, source clustering, photo-z into the forward model. Then stress-test peaks, the ℓ₁‑norm, the compressor.

Not specific to lensing: a trained denoiser in a proximal scheme fits any linear inverse problem with known noise. Galaxy clustering; 21‑cm, where foreground removal sits where mass mapping sits here.

the honest opening — everything in the thesis is simulations, on purpose, because isolating one analysis step needs the ground truth. Taking it to data is the next thing
▲ if you name only one — connect the two halves. The mass-mapping side and the statistics side were built in parallel and never run end to end. The thesis calls that the most direct continuation
the learned-map caveat, if pressed — a network trained at a fiducial cosmology cannot bias the inference when it is applied to the sims and the data alike. What it can do is lose constraining power away from its training point
the one that excites you — BNT-aware reconstruction. It would recover at the map level the cross-bin information that per-channel statistics cannot reach, which is exactly the puzzle Part 4 ended on
if asked about other probes — the thesis proposes the method transferring, not a cross-correlation analysis. Galaxy clustering first, 21-cm as the interesting one, because foreground removal occupies the same slot in the chain that mass mapping does here

source material — not part of the talk

Lensing pedagogy and mass-mapping methods

21 slides lifted unchanged from LAM 2026, for Act 0 and the Ch2 backup library.

$\,$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$

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.

From galaxies to mass maps

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).

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}$$

Trustworthy non-Gaussian inference for weak-lensing cosmology

The Universe evolved from quantum fluctuations to galaxies over 13.8 billion years

ΛCDM describes the Universe with six free parameters, plus a set of assumptions

Independent probes test the model through both geometry and growth

Foreground mass deflects the light of background galaxies and distorts their images

Euclid: Mapping the Dark Universe

Real weak-lensing cosmological analyses are extremely complicated

The analysis chain from galaxy shapes to parameters

The analysis chain from galaxy shapes to parameters, and the four questions of this thesis

The next two parts live in one step of that chain

Part 1: Impact of weak-lensing mass-mapping algorithms on cosmology inference

We measure the shear, but theory predicts the convergence

The convergence is a projection of the matter density, sensitive to geometry and growth

Shear and convergence are both second derivatives of the lensing potential

Relation between $\kappa$ and $\gamma$

In practice, mass mapping is an ill-posed inverse problem

Kaiser–Squires inversion is exact for complete, noiseless data, but amplifies the noise

Every mass-mapping method minimises a data term plus a regulariser that encodes the prior

MCALens models κ as a Gaussian plus a sparse component

Different priors give different maps. Does the difference affect the cosmological results?

§1 The simulations, statistic and likelihood are fixed, and only the mass-mapping method varies

§1 MCALens recovers small-scale structure that Kaiser–Squires smooths out

§1 MCALens improves the figure of merit by 157 % over Kaiser–Squires, by making more information in the data available.

§1 The gain comes from scales below 8′, which only MCALens recovers

Part 2

A plug-and-play approach with fast uncertainty quantification for weak lensing mass mapping

What about Deep Learning for Mass Mapping?

Deep Learning for Mass Mapping?

Plug-and-Play Mass Mapping

How do we estimate uncertainties?

§2 PnPMass is as accurate as the fine‑tuned networks, with the smallest calibrated error bars

Performance

Conclusions

LAM originals, Part 1

Mass mapping methods:

Inference pipeline

So does the choice of the mass mapping algorithm matter for the final constraints?

Inference results

Wavelet multi-scale peak counts

Where does this improvement come from?

Restyled Acts 1–2

Is mass mapping preprocessing, or does the choice of reconstruction change the cosmology we infer?

Can a reconstruction be flexible, fast, accurate and honest about its own uncertainty at the same time?

§2 Every existing reconstruction gives up at least one of the four

Wiener, MCALens

DeepMass, MMGAN

all four at once

§2 PnPMass: one denoiser, trained once, inside a fixed‑point iteration

§2 As accurate as a network retrained for the observation — with the smallest calibrated error bars of any method tested

Parts 3 and 4 move from the map to the summary statistic

The map has to be compressed into a summary statistic before inference

The two‑point function is the standard statistic, and is complete for a Gaussian field

Traditional approach: the 2-point route

Two very different fields can have the same power spectrum

Forecasts show that higher-order statistics constrain better than the power spectrum

Indeed, HOS can tighten constraints significantly

Peak funtions

Higher Order Statistics:Peak Counts

Wavelet peaks

Wavelet peak funtions

Starlet $\ell_1$-norm

Inference with HOS

Noise & Covariance

From data to contours

Summary statistics

Key concept: wavelets

Wavelet peak counts and the $\ell_1$-norm are one-point statistics of the starlet coefficients

The summary statistic is compared with theoretical predictions to infer the cosmological parameters, in a Bayesian framework

Classical inference needs an explicit likelihood, usually Gaussian in the data vector

Key concept: generative modelling

A normalizing flow maps a simple distribution to a complex one through an invertible transform

Simulation‑based inference learns the posterior directly from simulated data

Part 3

The joint wavelet ℓ1‑norm matches neural compression for tomographic weak‑lensing inference

How much of the map's information can a summary statistic keep?

A neural compressor trained to be information-optimal

§3 Both summaries are compared on the same maps, through the same flow, both calibrated

§3 The CNN beats the ℓ1-norm by 36% in the figure of merit

§3 But could we do better than that? Weak lensing tomography

§3 The cross-bin information can be reached through cross maps or through a joint statistic

The joint wavelet ℓ₁‑norm matches neural compression for tomographic weak‑lensing inference

The late-time S₈ and the local H₀ disagree with the CMB predictions

Learned vs analytical,
and can we trust it?