ELUCID-DESI I: A Parallel MPI Implementation of the Initial Condition Solver for Large-Scale Reconstruction Simulations

Feng Shi; Hong-Ming Zhu; Huiyuan Wang; Jiaxin Han; Junde Li; Qingyang Li; Wensheng Hong; Xiaohu Yang; Yizhou Gu; Youcai Zhang

REVIEW 5 major objections 5 minor 1 cited by

ELUCID-DESI I: A Parallel MPI Implementation of the Initial Condition Solver for Large-Scale Reconstruction Simulations

T0 review · 5 major / 5 minor · reviewed 2026-08-03 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read An MPI-parallelized Hamiltonian Monte Carlo reconstruction code with a fast initial-guess module makes 8192^3-particle cosmic initial-condition reconstructions feasible on current HPC hardware.

desk verdict A credible MPI port of the ELUCID HMC reconstruction code, but the burn-in savings metric is internally inconsistent and the 8192^3 cost projections rest on an unmeasured extrapolation. read the letter →

arxiv 2602.06463 v2 pith:RE74DFMC submitted 2026-02-06 astro-ph.GA

Wensheng Hong , Xiaohu Yang , Junde Li , Huiyuan Wang , Zhao Chen , Hong-Ming Zhu , Qingyang Li , Yizhou Gu

show 5 more authors

Youcai Zhang Feng Shi Jiaxin Han Yu Yu Zhongxu Zhai

This is my paper · ORCID

classification astro-ph.GA

keywords cosmicinitialconditionsHamiltonianMonteCarloMPIdomaindecompositionFastPMburn-inreductiondensityfieldreconstructionELUCID-DESI

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper claims that Bayesian reconstruction of the cosmic initial density field, previously limited by single-node memory and long burn-in phases, can be scaled to survey-sized volumes through two engineering advances. First, a distributed-memory MPI implementation with domain decomposition and a fast particle-mesh solver shows nearly linear scaling in both particle count and CPU core count, with memory per process following the ideal line. Second, a new initial-guess module starts the Markov chain from a fast approximate linear inversion of the observed density field, cutting burn-in by 53 steps for 256^3 particles and 106 steps for 1024^3 particles, saving thousands of core-hours. The authors argue that these gains together make an 8192^3-particle reconstruction schedulable on about 512 nodes, a scale previously out of reach.

What carries the argument

The load-bearing objects are: the Hamiltonian Monte Carlo sampler, whose potential energy is the negative log-posterior of the initial density field; the FastPM particle-mesh solver, whose modified drift and kick factors reproduce linear growth exactly and serve as the forward model; the three-dimensional domain decomposition that distributes particles and density grids across MPI processes; and the initial-guess module, which performs a fast approximate inversion in Fourier space, dividing by a transfer function built from emulator power spectra. A scale-dependent transfer function calibrates the fast solver against a high-accuracy reference so the likelihood remains unbiased, while the gue

What would settle it

Run a 2048^3 reconstruction with both random and guess initialization and compare the measured CPU hours to the extrapolated values: if per-step wall time grows superlinearly with particle count, or if the burn-in reduction is far below the projected roughly 180 steps, the feasibility estimates for 8192^3 would need to be revised downward.

Watch

Extended reading notes

Core claim

The central claim is that the two bottlenecks of Hamiltonian Monte Carlo initial-condition reconstruction, the single-node memory ceiling and the expensive burn-in period, can both be removed at scale. The MPI code distributes particles and grids across nodes via three-dimensional domain decomposition, and although communication overhead makes each step about eight times slower than the shared-memory version, wall-clock time still grows almost linearly with particle number and core count for large problems. The initial-guess module applies a transfer-function-weighted linear inversion of the observed final density field to produce a chain start that already lies near the high-probability reg

Load-bearing premise

The load-bearing premise is that the nearly linear scaling and burn-in savings measured at 256^3 to 1024^3 particles extrapolate linearly to 8192^3 particles with 32,768 MPI processes; the paper contains no runtime measurement beyond 1024^3.

Editorial extensions

If this is right

Bayesian initial-condition reconstruction at 8192^3 particles becomes schedulable on about 512 nodes with roughly 400 GB of memory per node.
The guess module's savings grow with particle count, so the relative cost of burn-in shrinks as survey volumes grow.
Removing the single-node memory ceiling allows the same HMC formalism to be applied to volumes about 64 times larger than earlier reconstructions.
Reconstruction fidelity is maintained: phase correlations match the established characteristic scale relation and reconstructed power spectra agree within about 5 percent for k > 0.05 h/Mpc.
For a 4096^3 DESI-scale reconstruction, the projected total cost is about 5 million core-hours, with the guess module saving roughly 2.3 million core-hours, about 45 percent of the total.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

If burn-in savings continue to grow with volume, the dominant cost shifts from initialization to the post-burn-in sampling phase, so future gains would come from faster forward solvers or GPU offloading rather than better starting points.
The measured eightfold per-step MPI slowdown suggests that a hybrid strategy, shared-memory parallelism within a node and MPI across nodes, could recover single-node efficiency without losing the memory advantage.
A direct 2048^3 burn-in test would be the natural check on the extrapolation; without it, the 37-million-core-hour savings figure for 8192^3 remains a projection, not a measurement.
The guess module could likely be pushed further with a nonlinear or learned inversion; if the starting chi-square drops further, convergence could become even faster than the linear-inversion estimate suggests.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

The MPI port is real, and the scaling measurements at 256^3–1024^3 look plausible. That is the core contribution: a distributed-memory version of the ELUCID HMC framework with domain decomposition, clean memory scaling, and honest accounting of the ~8x per-step slowdown relative to OpenMP. The guess module—a fast linear inversion as an HMC warm start—is a sensible idea, and the validation against Gadget-4, Kun, and Quijote is the right way to show the reconstruction still works. I believe the engineering is useful and the accuracy checks are solid.

The soft spots are in the performance claims. The burn-in reduction numbers in Table 1 rely on a definition that does not hold together. They define convergence as chi^2 dropping below 1% of its starting value. For the random-start chain, the threshold at 256^3 is ~0.0054, just above the plateau. For the guess-start chain, the starting chi^2 is much lower (the head-start numbers imply it equals the random chain's value after 27–90 steps), so 1% of that lower value is far below the plateau. By the stated definition, the guess chain would be judged to have converged much later, if ever, not 53–106 steps sooner. Either the definition is misstated or the savings are not robust. Since these savings are the basis for the 8192^3 feasibility estimate, this is a load-bearing flaw.

Second, the 720 steps and ~37 million CPU hours for 8192^3 are extrapolated from runs that stop at 1024^3 (with one 2048^3 point for the guess advantage). No measurement exists beyond 1024^3, and the abstract/full-text mismatch (4096^3 vs 8192^3 as the planned run) does not inspire confidence. The scaling plots also normalise the 1024^3 curve by dividing by 4 to force it onto the ideal line, and the particle-number scaling curves are admitted to be skewed by load imbalance at the smallest size. The quoted 'nearly linear' scaling is therefore more forgiving than the figures strictly show.

Finally, no code or data release for a methods paper is a letdown; 'available on reasonable request' is not reproducibility.

Who this is for: people working on Bayesian initial-condition reconstruction or planning large constrained simulations. It deserves a serious referee, but the burn-in definition must be corrected and the extrapolation re-derived before the headline cost savings can be trusted. If the authors redo the comparison using a fixed convergence threshold or demonstrate both chains reach the same posterior, the paper becomes solid.

Referee Report

5 major / 5 minor

Summary. The paper presents an MPI-parallelized implementation of the ELUCID Hamiltonian Monte Carlo (HMCMC) initial-condition reconstruction framework, built on the FastPM particle-mesh solver with 3D domain decomposition, together with an initial-guess module that generates an approximate linear inversion of the observed density field to shorten burn-in. The authors report strong- and weak-scaling tests for 256^3, 512^3, and 1024^3 particle boxes using up to 2048 MPI processes, memory-scaling behavior, reconstruction-accuracy comparisons against Gadget-4, Kun, and Quijote simulations, and burn-in savings for 256^3 and 1024^3 runs. From these measurements, they extrapolate CPU-hour savings for 2048^3 and 8192^3 runs and argue that survey-scale ELUCID-DESI reconstructions are computationally feasible.

Significance. If the central claims hold, this is a useful contribution to Bayesian initial-condition reconstruction: it removes the single-node shared-memory bottleneck, demonstrates nearly linear scaling in the tested range, and proposes a practical way to shorten the expensive burn-in phase. The accuracy validation against three independent simulation suites is a genuine strength, as is the explicit treatment of the transfer-function calibration and the FastPM integration. The measured results at 256^3 and 1024^3 are valuable even if the larger extrapolations are not accepted. However, the headline quantitative claims—especially the CPU-hour savings and feasibility of 8192^3 runs—depend on a convergence definition that is internally inconsistent and on extrapolations that are not supported by data in the manuscript. The paper is therefore technically promising but needs substantial revision before the quantitative conclusions can be relied upon.

major comments (5)

[Abstract vs. §3.5, Table 1] The abstract and the body report different target scales and different projected costs. The abstract states that the planned ELUCID-DESI reconstruction uses 4096^3 particles, requiring about 800 HMCMC steps (~5 million core hours), with about 360 steps (~2.3 million core hours) saved by the guess module. The full text, however, repeatedly refers to an 8192^3 planned run (§3.3.1, §3.5, Table 1) and estimates about 720 steps and ~37,000,000 core hours saved. These are materially different claims. The authors must reconcile the abstract with the body and specify which target is actually planned.
[§3.5, Table 1 — burn-in definition] The burn-in period is defined as the number of iterations for chi2_omega to drop below 1% of its starting value for the first time. This definition is applied to both chains, but the two chains start from very different chi2 values. For the 256^3 random-start chain the initial value is chi2 ≈ 0.54 (§3.4.1), so the threshold is ≈ 0.0054; the guess-start chain begins at a much lower value, equivalent to that reached by the random chain only after 27 steps (Table 1). If 1% of the guess chain's own starting value were used, the threshold would typically lie below the converged plateau (≈ 0.005 at 256^3), so the guess-start chain would not count as converged at all under the stated rule. The reported savings of 53 and 106 steps are therefore not a valid comparison under the stated definition. Please define burn-in relative to a common absolute criterion and report the starting chi2 and conver
[§3.5, Table 1 — extrapolated CPU-hour savings] The rows for 2048^3 and 8192^3 are projections, not measurements. No wall-clock timings are reported at 2048^3 in the scaling sections, and the 8192^3 estimate assumes that per-step runtime scales linearly with particle number and that the head-start in steps grows linearly with box size. The implicit per-step CPU-hour costs in Table 1 are not internally consistent: from 1024^3 to 2048^3 the particle number increases by a factor of 8, but the per-step CPU hours increase by a factor of about 11.4, whereas from 2048^3 to 8192^3 the particle number increases by a factor of 64 and the per-step cost increases by exactly 64. The hidden assumptions about core counts and node layouts need to be stated, and the 37 million CPU-hour savings should be presented as an idealized extrapolation rather than as a measured or nearly measured result.
[§3.2.1, Fig. 2 left — normalization of 1024^3] The 1024^3 processor-scaling curve is placed on the common ideal line by first normalizing to its runtime at 256 processes and then dividing by a factor of 4, assuming ideal scaling between 64 and 256 processes. No measurement at 64 processes exists for this problem size. This anchors the 1024^3 curve at the ideal line by construction. The slope from 256 to 2048 processes is measured, but the vertical offset is assumed. The visual statement that the 1024^3 curve deviates the least from ideal scaling is therefore partly an artifact of the normalization. Please show unnormalized wall-clock times or compute parallel efficiencies from measured anchors at a common process count.
[§3.2.2, Fig. 2 right — load-imbalance normalization] The right panel of Fig. 2 is used to support the claim of near-linear scaling with particle number. The text acknowledges that the apparent flattening of the slopes with increasing Nproc is caused by load imbalance at the 256^3 normalization point, which makes the normalized runtimes of larger problems look artificially small. This means the plot, as presented, is not a clean measurement of particle-number scaling. The authors should plot unnormalized runtimes or apply an explicit load-imbalance correction and report parallel efficiency values, so that the linear-scaling claim can be assessed directly from the data rather than through a known artifact.

minor comments (5)

[§2.4.1] The subsection heading 'Motivation and' appears incomplete; the missing part should be supplied.
[Eq. (28)] Equation (28) uses the notation rho_p in the inversion but only rho_obs was defined earlier. Please define all symbols consistently.
[Introduction] The citation 'Sousa et al. 2007' in the context of peculiar-velocity field reconstruction appears to be a mismatched reference; the listed source is not a velocity-field paper. Please check and replace it.
[Fig. 4 and Table 1] The 2048^3 case appears in Fig. 4 and Table 1, but the setup of that run (Nproc, node layout, solver parameters) is not described in Section 3.5. Please provide the run configuration or clearly mark the 2048^3 row as an extrapolation.
[Figure 5 caption] The caption lists 'Bottom left' and then 'Middle right' for the two bottom-row panels, which is confusing. The panel labels should be corrected.

Circularity Check

0 steps flagged · score 2.0 of 10

No significant circularity: core scalability and accuracy claims are externally benchmarked; the burn-in-savings metric is chain-dependent but not an input-output identity.

full rationale

The paper's central MPI scaling and reconstruction-accuracy claims are measured against wall-clock benchmarks and independent simulations (Gadget-4, Kun, Quijote), not derived from the model's own outputs. The transfer function T(k) is calibrated from simulation pairs and injected into the forward model as an input; it is not used as the success metric. The guess module is an approximate linear inversion (Eqs. 27-28) whose starting field is then evolved and scored by the same chi-squared likelihood; the reported acceleration is an empirical chain comparison, not an algebraic consequence of the guess formula. The one non-circular but load-bearing concern is in Section 3.5/Table 1: burn-in is defined as the first time chi-squared drops below 1% of its own starting value, making the comparison depend on each chain's initial point and making the extrapolated 2048^3/8192^3 CPU-hour savings fragile. That is an internal-consistency/robustness problem, not a case where the prediction reduces to the input by construction. Self-citations to Wang et al. (2013, 2014) and the CSST emulator are background methodology or calibration inputs, and the accuracy validation uses independent external simulations, so they do not create circularity.

Assumptions & free parameters 3 free parameters · 5 assumptions · 0 invented entities

The central performance claims rest on the Gaussian prior, a calibrated PM forward model, a self-defined convergence threshold, and untested linear extrapolation to 8192^3. No invented physical entities are introduced.

free parameters (3)

R_s (Gaussian smoothing scale) = 3.0 h^-1 Mpc for 256^3; 5.0 h^-1 Mpc for 1024^3
Chosen per run; controls the scale down to which phase information is recovered and affects chi2 and power-spectrum comparisons.
R_pm (PM force resolution) = 0.75 for 256^3; 2.0 for 1024^3
Chosen per run; sets the force resolution of the particle-mesh solver and affects reconstruction accuracy.
N_pmstep (number of PM time steps) = 10
Chosen as a default for all tests; directly sets the cost and accuracy of each forward model evaluation.

assumptions (5)

domain assumption Initial density field is a Gaussian random field with linear power spectrum P_lin(k).
Standard ΛCDM assumption used in Eq. (4); if it fails, the posterior is misspecified.
domain assumption FastPM/PM evolution plus transfer function T(k) adequately approximates a full N-body final density field.
T(k) is calibrated from paired simulations (Eq. 1-2); the small-box/emulator calibration in Appendix B assumes transferability to the full reconstruction volume.
ad hoc to paper Convergence is defined as chi2_omega dropping below 1% of its starting value for the first time.
Operational definition used to quantify burn-in savings; it is not a formal MCMC convergence diagnostic and can make speedup numbers look larger than they are.
ad hoc to paper Scaling and burn-in savings extrapolate linearly from 256^3-1024^3 to 2048^3 and 8192^3.
Table 1 linearly extrapolates steps and core-hours to scales that were never run; this is the load-bearing assumption for the flagship feasibility claim.
domain assumption Calibration data from a small-box simulation and the CSST emulator are valid for the full volume.
Appendix B assumes the transfer function and power-spectrum ratio are not volume-dependent; if false, guess quality and likelihood correction degrade.

how reviews work

0 comments

Cite this review

Pith. "Pith review of ELUCID-DESI I: A Parallel MPI Implementation of the Initial Condition Solver for Large-Scale Reconstruction Simulations." pith.science (2026). https://pith.science/paper/RE74DFMC

@misc{pith2026260206463,
  author       = {Pith},
  title        = {Pith review of: ELUCID-DESI I: A Parallel MPI Implementation of the Initial Condition Solver for Large-Scale Reconstruction Simulations},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RE74DFMC}},
  note         = {Machine review of arXiv:2602.06463}
}

abstract

We present a highly scalable, MPI-parallelized framework for reconstructing the initial cosmic density field, designed to meet the computational demands of next-generation cosmological simulations, particularly the upcoming ELUCID-DESI simulation based on DESI BGS data. Building upon the Hamiltonian Monte Carlo approach and the FastPM solver, our code employs domain decomposition to efficiently distribute memory between nodes. Although communication overhead increases the per-step runtime of the MPI version by roughly a factor of eight relative to the shared-memory implementation, our scaling tests-spanning different particle numbers, core counts, and node layouts-show nearly linear scaling with respect to both the number of particles and the number of CPU cores. Furthermore, to significantly reduce computational costs during the initial burn-in phase, we introduce a novel ``guess'' module that rapidly generates a high-quality initial density field. The results of the simulation test confirm substantial efficiency gains: for $256^3$ particles, 53 steps ($\sim$ 54 core hours) are saved, accelerating convergence by a factor of $\sim$ 18; for $1024^3$, 106 steps ($\sim$7500 core hours), achieving a speedup factor of $\sim$ 3. The total core hour gain grows with the number of particles, rendering large-volume reconstructions computationally practical for upcoming surveys, including our planned ELUCID-DESI reconstruction simulation with $4096^3$ particles. We estimate that achieving convergence for this scale (targeting DESI-BGS data) requires about 800 HMCMC steps ($\sim$ 5 million core hours). Our initial guess module will save approximately 360 steps ($\sim$2.3 million core hours), reducing the total computational time by about 45\%.

Figures

Figures reproduced from arXiv: 2602.06463 by the authors.

**Figure 1.** Schematic overview of the reconstruction pipeline. The workflow consists of two main stages. Left panel: Pre-processing stage, which prepares the necessary inputs: (a) reconstruction of the three-dimensional observed density field from survey data, and (b) calibration of the transfer function 𝑇 (𝑘) to correct for inaccuracies in the fast particle-mesh forward model. Right panel: The core HMCMC reconstruction loop, w… view at source ↗

**Figure 2.** Parallel scaling performance. (Left) Processor-number scaling for fixed particle numbers (2563 , 5123 , 10243 particles). Runtime per iteration is normalized to show parallel efficiency relative to ideal linear scaling (dashed line). Larger problems scale better. (Right) Particle-number scaling at fixed MPI process counts. Runtime scales nearly linearly with total particle count (dashed line). The flatter slopes for… view at source ↗

**Figure 3.** Memory usage per MPI process. (Left) Processor-number scaling: Memory per process for fixed reconstruction particle numbers (2563 , 5123 , 10243 ), normalized to show efficient reduction relative to ideal scaling (dashed line). (Right) Particle-number scaling: Memory per process as a function of total particle count, with the number of MPI processes held constant for each curve. All curves exhibit highly linear scal… view at source ↗

Figures from the paper (4 more)

**Figure 4.** Figure 4: 𝜒 2 𝜔 values as a function of the number of steps N. Results are shown for chains started from a random field (solid curves) and from the guess module (dashed curves) for three problem sizes (2563 , 10243 , 20483 ). In all cases, the guess module delivers a better init…

**Figure 5.** Figure 5: Visual and quantitative assessment of a 2563 reconstruction. Top left: True initial density field at 𝑧ini. Top middle: Reconstructed initial density field. Top right: Relative error map, showing no coherent large-scale bias. Bottom left: Final density field at 𝑧 = 0 ev…

**Figure 6.** Figure 6: Final density fields at 𝑧 = 0 for the 10243 test case. Left: Field obtained by evolving the reconstructed initial conditions with Gadget-4. Middle: Original input field from the Kun simulation. Right: The relative error map. The overall cosmic web structure is faithful…

**Figure 8.** Figure 8: Density-density scatter plot between the original and reconstructed 𝑧 = 0 density fields (Gaussian smoothed with 𝑅 = 4.8 ℎ −1Mpc). Contours enclose 68%, 95%, and 99% of grid cells. The tight correlation around the one-to-one line (solid black) indicates unbiased recons…

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

Multi-tracer mass bias in matched cosmic voids from SDSS DR7 and the ELUCID constrained simulation
astro-ph.CO 2026-08 conditional novelty 6.0 of 10

Galaxy and subhalo mass fractions relative to dark matter decrease towards the centres of matched cosmic voids, with galaxy-to-subhalo ratios limited by the scarcity of massive subhaloes.

Reference graph

Works this paper leans on

1 extracted references · cited by 1 Pith paper

[1]

mini-simulation

Bertschinger E., 1987, ApJ, 323, L103 Blanton M. R., et al., 2005, AJ, 129, 2562 BlasD.,LesgourguesJ.,TramT.,2011,JournalofCosmologyandAstropar- ticle Physics, 2011, 034–034 CSST Collaboration et al., 2025, arXiv e-prints, p. arXiv:2507.04618 Chen Z., Yu Y., 2025a, Extending CSST Emulator to post-DESI era (arXiv:2510.09503),https://arxiv.org/abs/2510.0950...

arXiv 1987

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.

[1] [1]

mini-simulation

Bertschinger E., 1987, ApJ, 323, L103 Blanton M. R., et al., 2005, AJ, 129, 2562 BlasD.,LesgourguesJ.,TramT.,2011,JournalofCosmologyandAstropar- ticle Physics, 2011, 034–034 CSST Collaboration et al., 2025, arXiv e-prints, p. arXiv:2507.04618 Chen Z., Yu Y., 2025a, Extending CSST Emulator to post-DESI era (arXiv:2510.09503),https://arxiv.org/abs/2510.0950...

arXiv 1987