HAMLET: HMC for Cosmological Field Reconstruction

• HAMLET is a Bayesian algorithm that jointly reconstructs the linear large-scale density and 3D velocity fields from raw survey observations.
• It integrates true distances as latent variables to correct lognormal biases, employing Hamiltonian Monte Carlo and GPU acceleration for drastic speed improvements.
• Validation on Cosmicflows-like mocks demonstrates high correlation in reconstructed fields and superior performance compared to traditional Gibbs-sampling methods.

HAMLET, the HAmiltonian Monte carlo reconstruction of the Local EnvironmenT, is a Bayesian forward-modeling algorithm for reconstructing the linear large scale density field $\delta(\mathbf r)$ and its associated 3D velocity field $\mathbf v(\mathbf r)$ from peculiar-velocity surveys such as Cosmicflows. Its defining feature is a joint posterior over Fourier density modes, true distances, and optional nonlinear-velocity dispersion, so that the lognormal distance–velocity transformation is treated inside the inference procedure rather than by an external pre-correction step. In the formulation introduced by Valade et al., HAMLET was tested on Cosmicflows mock catalogues with up to 30 000 data points and was reported to outperform earlier Gibbs-sampling reconstructions by two to four orders of magnitude in CPU time, with the gain attributed to both Hamiltonian Monte Carlo and GPU execution (Valade et al., 2022).

1. Reconstruction target and physical setting

HAMLET addresses the inference problem defined by peculiar-velocity cosmography: given measured distance moduli $\mu_i$ and redshifts $z_i$, infer the underlying large-scale density and velocity fields. In linear theory, the connection between density and velocity is expressed in Fourier space as

$$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$

so reconstruction of the Fourier modes $\tilde\delta(\mathbf k)$ determines $\mathbf v(\mathbf r)$.

The observational difficulty is not merely sparsity and noise. Distance moduli are modeled with Gaussian errors $\sigma_{\mu,i}$, but converting $\mu \rightarrow d$ and then to line-of-sight peculiar velocity introduces a lognormal bias, often conflated with Malmquist bias. HAMLET is designed to “undo” that effect by treating the true distances as latent variables in the posterior, instead of applying a separate correction to derived peculiar velocities. This places the distance transformation, the velocity prediction, and the cosmological prior inside a single inferential object (Valade et al., 2022).

A plausible implication is that HAMLET is best understood not as a post-processing method for catalogued peculiar velocities, but as a forward model from raw survey observables to latent cosmic fields. That distinction matters because the dominant non-Gaussianity enters through the distance transformation rather than through the prior on the linear density field.

2. Bayesian posterior, likelihood, and priors

The unknowns are the set of true distances $D=\{d_i\}$, the Fourier modes $\mathbf v(\mathbf r)$0, and optionally a nonlinear-velocity dispersion $\mathbf v(\mathbf r)$1. With data $\mathbf v(\mathbf r)$2, the posterior is

$\mathbf v(\mathbf r)$3

The likelihood factorizes over objects and over distance and radial-velocity components,

$\mathbf v(\mathbf r)$4

For distances, HAMLET uses

$\mathbf v(\mathbf r)$5

with $\mathbf v(\mathbf r)$6.

For radial peculiar velocities, the observed quantity is defined as

$\mathbf v(\mathbf r)$7

while the model prediction is

$\mathbf v(\mathbf r)$8

The corresponding likelihood is

$\mathbf v(\mathbf r)$9

with

$\mu_i$0

The prior factorizes as $\mu_i$1. The density prior is a Gaussian random field,

$\mu_i$2

where $\mu_i$3 is the linear matter power spectrum. The distance prior approximates the survey selection function through the histogram of observed redshift-distances, and $\mu_i$4 is taken to be weak, for example uniform or broad Gaussian (Valade et al., 2022).

This posterior architecture makes the role of the data model explicit. The density field remains Gaussian at the prior level; the nonlinearity enters through the observational map from modulus and redshift to distance and radial velocity.

3. Hamiltonian Monte Carlo formulation

HAMLET bundles the latent variables into a parameter vector

$\mu_i$5

The negative log-posterior defines the potential energy,

$\mu_i$6

Introducing conjugate momenta $\mu_i$7 of matching dimension and a tunable mass matrix $\mu_i$8, HAMLET defines the Hamiltonian

$\mu_i$9

The mass matrix is chosen, or periodically adapted, to approximate the covariance of $z_i$0. In the prior-only limit, the diagonal blocks can be set using the prior variances, including $z_i$1 for each $z_i$2 mode and an approximate variance $z_i$3 for each $z_i$4.

Sampling proceeds with the standard leapfrog integrator applied to Hamilton’s equations,

$z_i$5

Leapfrog is used because it guarantees reversibility and volume conservation. The integration error in $z_i$6 is $z_i$7 per step and does not accumulate with the number of steps. The trajectory length is $z_i$8, and the step size $z_i$9 is tuned to balance proposal quality against traversal efficiency; the implementation targets a Metropolis acceptance rate of approximately 0.65. Dual averaging can tune $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$0 during warm-up, and the No-U-Turn Sampler can optionally adapt $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$1 to avoid back-and-forth trajectories (Valade et al., 2022).

A common misconception is that HMC in this setting is simply a faster random-walk sampler. HAMLET’s formulation is more specific: the computational gain comes from using posterior gradients of a structured cosmological forward model, so the efficiency improvement is tied to geometry-aware proposals rather than to generic parallelization alone.

4. GPU execution and computational profile

HAMLET implements the expensive operations of the sampler through TensorFlow on GPUs. The documented GPU-resident components include FFT, the gradient of $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$2, and leapfrog updates. The gradient combines analytic derivatives of the log-likelihood and log-prior with an FFT-transpose for the $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$3 real-space transforms.

For a typical grouped CF3 problem with approximately $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$4 data points and a $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$5 Fourier-cell grid, HAMLET is reported to run a full HMC chain of approximately $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$6–$$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$7 steps in approximately 10 minutes on a single high-end GPU, compared with weeks on one CPU for the earlier Gibbs-sampler implementation. The corresponding acceleration is stated as approximately $$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$8–$$\tilde{\mathbf v}(\mathbf k)= -i\,H_0\,f(\Omega_0)\,\frac{\mathbf k}{k^2}\,\tilde\delta(\mathbf k),$$9. The abstract summarizes the same result as an improvement of two to four orders of magnitude in CPU time, attributed both to the higher efficiency of HMC and to GPU rather than CPU execution (Valade et al., 2022).

The performance claim is significant because the parameter space is not small. In the mock-catalogue tests, HAMLET is described as exploring approximately 50 times more dimensions than the prior Gibbs sampler while remaining faster in wall-clock time. This suggests that the principal scalability bottleneck shifts from Markov-chain inefficiency to memory and transform throughput.

The reported memory scaling is correspondingly favorable: the GPU-HMC framework is said to handle approximately $\tilde\delta(\mathbf k)$0 or $\tilde\delta(\mathbf k)$1 Fourier grids—about $\tilde\delta(\mathbf k)$2–$\tilde\delta(\mathbf k)$3 modes—at approximately 1–2 GB GPU memory. This is presented as sufficient for Cosmicflows-4, which is expected to contain approximately 5 times more velocities, together with higher-resolution modes.

5. Validation on Cosmicflows-like mock catalogues

The mock-catalogue campaign used nine mocks generated with $\tilde\delta(\mathbf k)$4, $\tilde\delta(\mathbf k)$5, or $\tilde\delta(\mathbf k)$6 galaxies inside $\tilde\delta(\mathbf k)$7, with distance-modulus errors set to $\tilde\delta(\mathbf k)$8 for $\tilde\delta(\mathbf k)$9. The underlying density and velocity fields came from a $\mathbf v(\mathbf r)$0 linear realization with Planck cosmology. Chains of approximately 200 warm-up steps followed by 500–1000 production steps produced unbiased estimates of $\mathbf v(\mathbf r)$1 and $\mathbf v(\mathbf r)$2, and convergence of the volume-averaged $\mathbf v(\mathbf r)$3 reached the percent level by approximately 100–200 steps (Valade et al., 2022).

Several reconstruction-quality diagnostics were reported. Mean-field slices $\mathbf v(\mathbf r)$4 reproduced the zero-contours of the target field out to approximately $\mathbf v(\mathbf r)$5–$\mathbf v(\mathbf r)$6 for CF3-like errors, and farther for $\mathbf v(\mathbf r)$7 or $\mathbf v(\mathbf r)$8. The local fractional uncertainty $\mathbf v(\mathbf r)$9 increased with radius, while smaller $\sigma_{\mu,i}$0 yielded tighter constraints. The radial profile of the bulk flow and monopole matched the target within the mock-boundary zone, $\sigma_{\mu,i}$1.

The radial degradation is quantifiable. Pearson correlations between the HAMLET mean field and the target decreased from approximately 0.9 at small radii to approximately 0.5–0.8 at $\sigma_{\mu,i}$2, with the correlations higher for velocities than for densities. This suggests that, within the linear-theory forward model adopted by HAMLET, peculiar-velocity data constrain low-order flow moments more robustly than small-scale density detail.

6. Comparative assessment, limitations, and uses

A subsequent comparative study tested HAMLET against BGc/WF and an “exact” Wiener filter on Cosmicflows-3-like mocks. In the nearby regime, $\sigma_{\mu,i}$3, the two practical methods were found to perform roughly equally well. HAMLET performed slightly better in the intermediate regime, $\sigma_{\mu,i}$4. The most substantial differences appeared in the distant regime, $\sigma_{\mu,i}$5, near the survey edge: HAMLET outperformed BGc/WF in terms of better and tighter correlations, but also yielded a somewhat biased reconstruction there, whereas such biases were reported to be absent from BGc/WF (2209.05846).

The distant-regime behavior is central to interpreting HAMLET’s output. The improvement in correlation does not imply uniformly unbiased reconstruction across radius. The testing study explicitly associates the far-edge degradation with a boundary effect, and this suggests that posterior richness alone does not remove the need for careful radial modeling of the survey boundary.

Within those limitations, the algorithm has two immediate applications. First, its acceleration makes possible larger reconstructions from upcoming Cosmicflows-4 data. Second, constrained realizations drawn from the HAMLET chain can be used as initial conditions for zoom-in N-body or hydrodynamic simulations of the Local Universe, and the abstract also notes the role of the method in setting constrained initial conditions for cosmological high resolution simulations (Valade et al., 2022).

The framework is also extensible. The parameter vector $\sigma_{\mu,i}$6 may be enlarged to include cosmological parameters such as $\sigma_{\mu,i}$7, $\sigma_{\mu,i}$8, and $\sigma_{\mu,i}$9, or nuisance parameters such as zero-points for different distance indicators. This suggests a broader role for HAMLET as a joint sampler over cosmography and calibration, rather than only a fixed-cosmology reconstruction engine.