Physics-Informed Machine Learning Workshop, FORTH, Heraklion, 16 September 2026

Machine Learning for Non-Gaussian Weak-Lensing Cosmology

Andreas Tersenov   FORTH, University of Crete, CEA Paris-Saclay

with Hubert Leterme, Jean-Luc Starck, Martin Kilbinger and Jalal Fadili

slides atandreastersenov.github.io/talks/PIML_FORTH_2026

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.

We measure a noisy, masked distortion field. The image we want is the projected mass

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
the image we want: 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)
the measurement: 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

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$, with $\hat{P}_1 = (k_1^2 - k_2^2)/k^2$, $\hat{P}_2 = 2k_1k_2/k^2$ in Fourier space
  • from shear to convergence: $\kappa = \hat{P}_1\gamma_1 + \hat{P}_2\gamma_2$, since $\hat{P}_1^2 + \hat{P}_2^2 = 1$
  • linear and exact for a complete, noiseless field; $\gamma = \mathbf{P}\kappa$ from here on

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

The analysis chain from galaxy shapes to parameters, and where the learning goes

question 1, the first half: the inverse problem

Learn only the prior. As accurate as end to end, any configuration, certified error bars?

question 2, the second half: the features

A learned encoder, or hand‑crafted features? Who extracts more about the parameters, and why?

Part 1

Learning the prior: an ill‑posed inverse problem, and what a learned prior buys


Tersenov, Baumont, Starck & Kilbinger 2025, A&A 698, A25
Leterme, Tersenov, Fadili & Starck 2026, A&A 710, A292
Leterme, Tersenov & Starck, in preparation

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

the measured shear field
measured shear
$\kappa = \;?$
a smooth, heavily denoised convergence map
a noise-dominated 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}$

the regulariser $\mathcal{R}$ encodes the prior on $\kappa$

1Kaiser–Squires no prior, smoothing only
2Wiener filtering an ℓ2 prior: a Gaussian field
3sparse recovery an ℓ1 prior: sparse in wavelets
4MCALens Gaussian + sparse
5deep learning a learned $\mathcal{R}$, from simulations

Kaiser & Squires 1993, Bobin et al. 2012, Lanusse et al. 2016, Starck et al. 2021, Jeffrey et al. 2020, Remy et al. 2023

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

simulations
cosmo‑SLICS
shear maps
Euclid‑like noise and mask
a heavily smoothed reconstruction
a noise-dominated reconstruction
a sharp reconstruction
mass maps
KS, iKS, MCALens
the only thing that varies
summary statistic
set of non‑Gaussian statistics
$\log \mathcal{L}(\theta)$
inference
posterior
Ωm, σ8, h, w0
the output we compare

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

§1 The sparse solver improves the map error by 4 %, and the constraining power by 157 %

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 reconstruction error by only 4 %, but the constraining power by 157 %
  • map error and constraining power are not the same objective: choose the reconstruction by the posterior it gives

Four ways to put a network in a linear inverse problem

end to endDeepMass, MMGANtrained on pairsdatanetworkimagenew mask or noiseretrain unrolledMRI, CTtrained throughnew mask or noisewithin training plug-and-playPnPMass, this talkdata stepdenoiser×8trained once, alonenew mask or noisenothing to retrain posterior samplingRemy et al. 2023data stepscore model×103samplestrained once, alonenew mask or noisenothing to retrain

Dashed boxes: what is trained together. Only the two on the right work for a new mask or noise covariance without retraining.

Plug-and-play: a learned denoiser replaces the proximal step

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{P}\kappa^n - \gamma \right) \right)$
  • The series converges to a fixed point $\hat{\kappa}$
  • choosing $\mathbf{B} = \mathbf{P}^T \Sigma^{-1}$ makes the method independent of the noise covariance and the mask

Venkatakrishnan, Bouman & Wohlberg 2013, Kamilov et al. 2023, convergence: Pesquet et al. 2021, Hurault et al. 2022

§1 A variant with more physics in: a Wiener filter takes the Gaussian part, the loop runs on the residual

the measured shear goes through a Wiener filter and gives the Gaussian component of the map; the ground truth beside it the Gaussian component's prediction is subtracted from the data; the residual goes through the plug-and-play loop with a denoiser trained on non-Gaussian residuals, and gives the residual estimate beside its ground truth

within 0.5 % of the end‑to‑end network, and converged in two iterations

Leterme, Tersenov, Fadili & Starck 2026

§1 A second network predicts the error, pixel by pixel, from the same simulated pairs

datadenoiser in the loop, ×8variance networksame simulated pairsthe map, κ̂one forward passits predicted error, σ̂one pass, no sampling
the objective: 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]$$

its minimiser is the posterior variance, pixel by pixel: one forward pass, no sampling

§1 Conformal calibration makes the error bar honest, on held‑out maps

schematic: one pixel across held‑out maps, at a visible miss rate 1. the network’s interval-20250 % of the truths outside 2. score, take the quantile0outsideinsideqhow far outside, sorted 3. widen every interval by q-20220 % outside: the target rate

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

for any network and any distribution, with n held‑out maps: distribution‑free, finite‑sample

Romano, Patterson & Candès 2019, Leterme et al. 2026

§1 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
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
a PnPMass convergence map beside its per-pixel calibrated uncertainty, on the same field
the map, and its calibrated σ

DeepMass and MMGAN are retrained per configuration; PnPMass is trained once, for any mask and noise level.

§1 Slice the sources by distance: six correlated channels, each measured with a sixth of the galaxies

channel 1z ∈ [0, 0.57]channel 2z ∈ [0.57, 0.80]channel 3z ∈ [0.80, 1.05]channel 4z ∈ [1.05, 1.36]channel 5z ∈ [1.36, 1.79]channel 6z ∈ [1.79, 2.60]the maps, per slice. the same structures in every channel, stronger with distancethe measurements. a sixth of the galaxies per channel, and many pixels with none

§1 Nulling: a fixed, invertible re‑mixing of the channels makes each one local in distance

before: channel k sees all the matter in front of itnearfarafter the nulling: each channel local in distance123456nearfarκ′ = B κlower‑triangular, invertiblethree channels at a timeweights from the geometry alonea difference of noisy maps

why: a structure’s apparent size mixes its physical size with its distance, so a multiscale analysis needs channels that are local in distance. The price: the noise adds up, and an error in one channel leaks into the next.

Bernardeau, Nishimichi & Taruya 2014; Barthelemy et al. 2022

§1 One multichannel denoiser in the same loop: the physics stays per channel, the correlation lives in the prior

per channel: six loops, six denoiserschannel kdata stepdenoiser kmap kiterate× 6, nothing sharedjoint: one loop, one denoiser, six channelssix channelsdata stepper channeldenoiser6 in, 6 outsix mapsiterateone problem, one prior over the six channels

$$\kappa^{n+1} = F_{\Theta}\!\left(\kappa^{n} - \tau\,\mathbf{B}\,(\mathbf{P}\kappa^{n} - \gamma)\right),\qquad \mathbf{P} = \mathrm{diag}(\mathbf{P}_0,\dots,\mathbf{P}_0),\qquad F_{\Theta} : \mathbb{R}^{6N^2} \to \mathbb{R}^{6N^2}$$

the same iteration: the operator, the masks and the noise levels enter only in the data step. New is the prior over the six channels together, and the theory follows it: convergence with masked data, and an error bound that fixes the step size.

§1 Joint reconstruction wins in every channel, and it is the only one that survives the foreground nulling

1. error per channel0.91.01.1the zero map123456per channeljointjoint 0.89, per channel 0.952. null the foregroundbefore: nestedafter: local in distancenearfardifferences of noisy maps3. after the nulling0.91.01.1the zero map1234561.48 on averagejointjoint 0.94, per channel 1.48

channel 1 of one field, the nearest slice

the measured shear of channel 1: noise and mask
measured
the per-channel reconstruction of channel 1: flat, nothing recovered
per channel
the joint reconstruction of channel 1: the peak and the structure below it
joint
the truth of channel 1: the same peak, faint
truth

Per channel sees nothing. Joint reads the peak from the other five channels, and it is in the truth.

Leterme, Tersenov & Starck, in preparation; κTNG, the nulling of Bernardeau, Nishimichi & Taruya 2014

Part 2

Learning the features, or not: hand‑crafted features against a learned encoder


Tersenov, Starck & Kilbinger, submitted to A&A

Two very different fields can have the same power spectrum

a large-scale-structure map, a Gaussian map, and their identical power spectra the same large-scale-structure map, then rebuilt from its Fourier phases only, and from its Fourier amplitudes only: the phases keep the filaments and haloes, the amplitudes keep nothing recognisable

The missing information is in the Fourier phases. To reach it we need statistics beyond second order: peak counts, wavelet statistics, the ℓ1‑norm, Minkowski functionals.

The $\ell_1$-norm is a one-point statistic of the wavelet coefficients, per scale and amplitude bin

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, one per angular scale
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$$

sum of absolute coefficients per band $j$ and amplitude bin $i$: every pixel counts

the feature vector
schematic l1-norm per signal-to-noise bin, one curve per starlet scale

one curve per scale, binned by amplitude (weighted histogram of the field in the wavelet domain)

The posterior is learned from simulator samples, likelihood‑free, and calibration‑tested

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

A learned encoder, trained jointly with the flow to maximise the mutual information \(I(t;\theta)\)

simulator
\(\theta\sim p(\theta)\)
→
example convergence-map patch
map \(x\)
→
encoder
\(f_\phi\)
→
summary
\(t\in\mathbb{R}^{10}\)
→
flow
\(q_\psi(\theta\mid t)\)
→
posterior
\(p(\theta\mid t)\)
\[\max_{\phi,\psi}\;\mathbb{E}_{p(\theta,x)}\!\left[\log q_\psi\big(\theta\mid f_\phi(x)\big)\right]\;\le\; I(t;\theta) - H(\theta)\]
parameter space  \(\theta\)

§2 The learned encoder is 36 % ahead of the hand-crafted features

matched-pipeline posterior for the auto-only l1-norm, with its figure of merit at 2448 the CNN posterior added: tighter, with a figure of merit of 3326 against the l1-norm's 2448
same maps, same flow, both calibrated
the gap
  • The learned encoder is 36 % ahead in constraining power
  • But the the ℓ1-norm is computed channel by channel → never sees the correlations between the maps

§2 The channels are correlated: each distance slice sees the matter in front of it

observer and the matter distribution the source galaxies as tomographic convergence maps, one per redshift bin
Credit: Justine Zeghal

§2 Two routes to the cross‑channel information: product channels, or a joint 2‑D histogram

route 1: product channels
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$

Strong only where both channels have structure. One map per pair, the same ℓ1‑norm.

route 2: a joint 2‑D histogram
channel j coefficient the joint plane of wavelet coefficients for a bin pair, with the two one-dimensional histograms along its edges the per-channel
ℓ1-norm only sees
the two marginal histograms channel 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 its pixels

§2 Read the channels jointly, and the hand‑crafted ℓ1‑norm matches the optimal encoder

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

Three open questions, on the steps where they live

coverage where it matters

The conformal guarantee is marginal, and the misses sit at the peaks. Conditional coverage at map level, at a cost a survey can pay?

certificates for large denoisers

Convergence needs a non‑expansive denoiser; for seven million parameters we verify it empirically. A cheaper certificate than a Jacobian bound?

learned features under a wrong simulator

The encoder's ceiling is the simulator's. When the physics is missing, which feature degrades gracefully, and can we tell without the truth?

Conclusions

  • Q1 ✓Learn the prior, keep the physics. Plug‑and‑play with a learned denoiser is as accurate as the end‑to‑end networks, its error bars are calibrated, and it is trained once, for any mask, noise level, or number of channels.
  • Q2 ✓Hand‑crafted features can match a learned encoder. Read the channels jointly, and the wavelet ℓ1‑norm ties the information‑optimal encoder, with no training.
  • so Learn only what the physics cannot supply, benchmark it, and calibrate it before you trust it.

Backup

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.

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

§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 The gain comes from scales below 8′, which only MCALens recovers

Kaiser–Squires 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
MCALens 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.

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)
Wavelets of different scales
  • 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
backup — not part of the talk

PnPMass

The denoiser and the transformer it is built from, the per-pixel uncertainties and their conformal calibration, the implementation, and the tomographic extension: its theory, the maps after the nulling, the paper’s figures.

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{P}^{\!\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.

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 tomographic loop: convergence with masked data, and an error bound that fixes the step size

Convergence, relaxed for masks

On the space of zero-mean maps (a centring layer at the output): the denoiser non-expansive with constant $L_\tau \le 1$, the data-step operator $\tilde{\mathbf{B}}\tilde{\mathbf{M}}\tilde{\mathbf{P}}$ symmetric positive semidefinite, and

$$0 < \tau L_\tau < 2\,\|\tilde{\mathbf{B}}\tilde{\mathbf{M}}\tilde{\mathbf{P}}\|^{-1}$$

give linear convergence to a unique fixed point at rate $L_\tau\rho_\tau$, with $\rho_\tau = \max(|1-\tau\lambda_{\min}|,\,|1-\tau\lambda_{\max}|)$. Masked pixels stay in: $\lambda_{\min} = 0$ here.

The error bound, and the step size

$$\mathbb{E}\|\hat\kappa - \kappa\|^2 \;\le\; \frac{1}{(1 - L_\tau\rho_\tau)^2}\;\mathbb{E}\|F_\Theta(\kappa + \mathbf{n}_\tau) - \kappa\|^2$$

The loop's error is the denoiser's own error on white noise of level $\tau$, times a factor that blows up as $\tau \to 0$. The safe choice is its first minimiser, $\tau = 2/\lambda_{\max}$: here $\lambda_{\max} = 7.66$, $\tau = 0.261$, the value the 2026 paper had picked by hand.

the normalised error of the joint reconstruction against the iteration number: from one at the start to about 0.89 after ten iterations, flat by twenty-four
24 iterations, against 8 in the single-map problem: a larger step size, a contraction rate closer to one

Leterme, Tersenov & Starck, in preparation, Propositions 1 and 2

After the nulling, the per‑channel reconstruction invents structure; the joint one stays quiet

channel 3z ∈ [0.80, 1.05]channel 4z ∈ [1.05, 1.36]channel 5z ∈ [1.36, 1.79]channel 6z ∈ [1.79, 2.60]per channelinvented structure,sign-reversed patchesjointquiet: close to thezero map,recovered barely

The same field as the main line, channels 3 to 6, one colour scale. The nulled channels carry little signal in this crop; the per‑channel maps fill them with foreground that the nulling should have removed.

Leterme, Tersenov & Starck, in preparation

The paper’s figures: error per channel for the denoiser alone, and inside the loop

normalised error per channel for the joint and the per-channel denoisers at noise level 0.2, before and after the nulling
denoising alone, white noise σ = 0.2 (about two galaxies per pixel)
normalised error per channel for the joint and the per-channel plug-and-play reconstructions after 24 iterations, before and after the nulling
inside the loop, 24 iterations

Blue joint, orange per channel; solid before the nulling, dashed after. Test set of the draft; the numbers are preliminary.

Leterme, Tersenov & Starck, in preparation, Fig. 4

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.

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
a convergence map convolved with a small wavelet gives a map of the small-scale structure; the same map convolved with a large wavelet gives a map of the large-scale structure

why wavelets for convergence maps

  • the LSS consists of localised structures (clusters, filaments, voids) with a characteristic size and position
  • scales can be analysed, or removed, separately

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 ℓ1‑norm, Minkowski functionals

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,...
generated images over a decade: a blurred grey 2014 face, sharper 2015 and 2016 faces, a photorealistic 2017 face, and a present-day generation of a fat ginger cat in a floral beret painting a portrait

Learn $\mathbb{P}_\theta$, then sample from it.

2014–2017 panels adapted from Shirvani, Image Generation Models: A Technical History, arXiv:2603.07455

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

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

§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 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), σ(σ8), σ(w0) and FoM3 across 9,000 mocks: medians 3045, 3371, 3326

The answer to “FoM3 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.

§2 Same maps, same flow, both calibrated: only the feature vector differs

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

Do baryons break HOS?

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

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.

optimistic pessimistic large, linear scales small, non-linear scales more → (schematic, not to scale) information beyond two-point baryonic feedback

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?

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

cosmoGRID kappa-maps 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)\)
JAX

§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

§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

§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
the CNN posterior in the standard and the nulled frame, the two contours lying on top of each other
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 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.

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 ℓ1‑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.