REVIEW 5 major objections 5 minor 40 references
Differentiable simulations for particle tracking in accelerators: analysis, benchmarking and optimization
T0 review · 5 major / 5 minor · reviewed 2026-08-06 · deepseek-v4-flash
Pith's one-line read Auto-differentiation beats finite differences for beamline gradients
desk verdict Useful first real benchmark of AD vs FD for accelerator tracking, with a clean complexity model; the headline wall-clock claim needs a parallel-FD CPU baseline before I'd trust the speedup factors. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is reverse-mode auto-differentiation (the adjoint method) applied to the directed acyclic graph formed by the elementary operations of a particle-tracking simulation. The paper's cost model expresses total gradient cost as $C_A/C_P = 1 + \kappa_{\mathrm{dp}} + \kappa_{\mathrm{sp}} + \kappa_{\mathrm{mp}} m$, where $m$ is the number of parameters and the coefficient $\kappa_{\mathrm{mp}}$ is small to negligible, in contrast with finite differences whose cost is $(am+b)C_P$. This makes the gradient cost effectively flat as parameters grow, and the paper validates the AD sensitivities against central finite differences before timing them.
What would settle it
Run reverse-mode AD gradient computation on a storage-ring lattice with thousands of turns and a few hundred parameters, and measure wall-clock time as a function of turns: if the time per gradient grows roughly linearly with the number of turns (or memory is exhausted) instead of staying near the cost of a few primal simulations, the central speedup claim would not carry over to lattice-scale problems.
Extended reading notes
Core claim
The central discovery is that reversing derivative propagation through the tracking simulation—reverse-mode auto-differentiation—makes the gradient cost nearly independent of the number of design parameters, while finite differences costs grow linearly with that number. The paper derives this asymmetry for the adjoint method, implements it without hand-deriving adjoint equations, and verifies it empirically: on CPU, AD-REV is about three times faster than finite differences on a single FODO cell and fifteen times faster on five cells; on a realistic 39-element, 58-parameter beamline it is an order of magnitude faster on CPU. Where GPU parallelism makes AD-REV slower on small problems, forward-mode AD is still faster than finite differences, so auto-differentiation wins on every tested configuration. A final optimization comparison shows gradient-based methods using AD gradients outperform Nelder-Mead both in speed of convergence and in the final objective value.
Load-bearing premise
The central speedup rests on the assumption that the computational graph of the tracking simulation can be held in memory during the backward pass; for deep lattices with many turns this memory cost could erase the advantage, a limitation the paper explicitly flags in Appendix B.
Editorial extensions
If this is right
- Beamline designers can swap finite-difference gradients for AD gradients inside standard optimizers and expect the per-iteration cost to stop growing linearly with the number of knobs, making mid- to high-dimensional design spaces tractable.
- The benchmarking protocol—cost model plus CPU/GPU timings normalized by primal cost—gives practitioners a way to predict which AD mode will be faster for a given lattice, particle count, and hardware.
- The optimization result indicates that gradient-based methods with AD gradients should become the default for transfer-line tuning, since they beat Nelder-Mead both in speed and in the quality of the final design.
- The advantage grows with the number of cells and parameters: on CPU the speedup over finite differences goes from roughly 3× on one FODO cell to 15× on five cells.
Reading between the lines
- The paper's Appendix B warns that reverse-mode AD must keep the whole computational graph in memory, so for storage-ring lattices with hundreds or thousands of turns the memory overhead could dominate; this implies forward-mode AD may be the better choice there, a step the paper does not take.
- The GPU results hint that a mixed strategy—forward-mode on small batches, reverse-mode on large—could recover performance across regimes, though the paper only reports the tradeoff.
- A direct testable extension: benchmark AD-REV on a realistic storage-ring lattice with turn-by-turn maps and measure whether wall-clock time grows linearly with turns (or memory exhausts), which would probe the Appendix B caution.
- From the paper's Eq. (11), one can extract a rule of thumb: reverse mode wins once the number of parameters is large compared with the primal's internal complexity; a dedicated experiment confirming this crossover would give designers a decision rule.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. This paper develops a JAX-based differentiable particle-tracking code and compares auto-differentiation (forward and reverse) with finite differences for gradient computation, both theoretically and empirically, on a space-charge FODO cell and a realistic CERN AWAKE transfer line. It reports CPU and GPU timings, validates AD gradients against central finite differences, and embeds reverse-mode AD gradients in SLSQP and L-BFGS-B optimizers, comparing them with Nelder-Mead. The central claims are that AD-REV is almost always faster than FD, that auto-differentiation is more efficient than finite differences in all tested cases, and that gradient-based optimization with AD gradients reaches a better minimum faster than gradient-free optimization.
Significance. If the empirical claims are robust, this is a useful first systematic benchmark of differentiable particle tracking, with a clearly derived complexity model and a welcome attempt to quantify the time-complexity distinction. The paper is candid about limitations, especially the memory cost of reverse-mode AD for deep lattices in Appendix B and the embarrassing parallelizability of finite differences in Section II.D. The code availability, the validation against finite differences, and the use of a realistic CERN beamline are concrete strengths. However, the headline conclusions are wall-clock claims, and the current evidence base lacks statistical rigor and an explicit treatment of the CPU finite-difference baseline, which makes the central speedup claims somewhat fragile.
major comments (5)
- [Section III.A.1, Section III.B, Figures 2-4 and 8] The central empirical claims are based on execution-time measurements reported without error bars, repeated runs, or any statement about averaging or run-to-run variability. For example, the speedup factors of 3x, 15x, and 10x in Figures 3 and 8 appear to be single point estimates. Since the conclusion that "AD-REV is almost always faster than FD" is a wall-clock claim, the manuscript should report statistics such as the median and spread over at least several independent runs, or justify why a single run is representative.
- [Section II.D, Figures 2a, 3, 8a] The finite-difference baseline on CPU is not described in terms of parallelization. Section II.D explicitly states that the perturbed primals of finite differences are embarrassingly parallel, and the reported CPU machine has 8 cores/16 threads. If the CPU FD measurements used a serial loop over the perturbed primals, the reported speedups may be substantially overstated relative to a reasonably optimized parallel FD implementation. The paper should state whether any CPU parallelization was used for FD and, if not, include a multithreaded or vectorized FD baseline, or justify why a serial FD is the appropriate comparison for the stated practical conclusion.
- [Section III.A.1, footnote 27] The manuscript reports that the AD gradients are "very similar" to those of JuTrack but states in footnote 27 that "the reason for the slight discrepancies is unknown." For a paper whose contribution is benchmarking and method comparison, an unexplained discrepancy with the published external reference weakens the validation claim. Please quantify the discrepancy (for example, relative error per parameter) and investigate likely sources, such as integration slicing, space-charge treatment, or finite-difference step size, or else soften the external-validation claim to a qualitative agreement with a stated tolerance.
- [Section III.B.1, Figure 9] The optimization comparison reports only one trajectory per optimizer with no repeated runs, seeds, or stated initial-simplex and stopping criteria. The conclusion that gradient-based methods "achieve their optima faster" and "reach a better minimum" is therefore a single-instance observation, not a robust comparison. Please add repeated optimizations from multiple starting points or random seeds, and report the distribution of final objectives and computation times, especially because Nelder-Mead and L-BFGS-B can be sensitive to initialization and numerical tolerances.
- [Conclusion and Appendix B] The conclusion states generally that "AD-REV is almost always faster than FD," but Appendix B itself notes that for lattices with hundreds to thousands of turns, AD-REV's memory cost for storing the computational graph may make it disadvantageous and that AD-FWD may be preferable. The tested cases are beamline-scale problems. The conclusions should be explicitly scoped to the beamline problems studied, and the abstract's reference to "beamlines and lattices" should be tempered, or additional evidence for lattice-scale applicability should be provided.
minor comments (5)
- [Appendix A.1, Figure 11, Eq. (A13)] There is an inconsistency in the reverse-mode AD example: for u3 = u1/u2 with u2 = 2, the partial derivative ∂u3/∂u1 is 1/u2 = 1/2, not u2. The accumulation step labeled "u8 · u2" and the vector entry "u2u8" in Eq. (A13) do not match the final analytical expression (1/2)(y^2 - x^2). Please correct the graph labels and the accumulation expressions.
- [Footnotes a and [25]] The code links are given as "can be found here" without showing a visible URL or repository identifier. Please include the actual repository URL or DOI in the text and in the bibliography so the code is accessible and the reproducibility claim is verifiable.
- [Figures 3, 4, and 8] The color bars in these figures show ratios of FD time to AD time, but the captions do not explicitly state that values greater than one mean AD is faster. Please state this in each caption to avoid ambiguity.
- [Figure 5] The x-axis label "FLOPS ×10^8" appears to denote the number of floating-point operations, not the rate in flops/s. Please relabel it, for example as "FLOPs (×10^8)", to be consistent with standard usage.
- [General] The phrase "almost always" in the conclusion is vague. Please replace it with a precise statement of the tested parameter ranges and the conditions under which FD was faster than AD-REV, such as small particle counts on GPU in Figures 2b and 4a.
Circularity Check
No significant circularity: the paper's gradient comparisons are validated against finite differences and an external reference, and its complexity model is derived analytically rather than fitted.
full rationale
The paper's central claims are empirical benchmarks comparing auto-differentiation (AD) with finite differences (FD), plus an analytical complexity model. The AD gradients are validated against central finite differences (Figs. 1 and 7) and against the published JuTrack values [18], which is an external benchmark rather than a self-citation or a fitted input. The complexity model in Eqs. (4)-(11) and Appendix A 2 is derived from first-principles operation counts (CP, CD, CC, κdp, κsp, κmp) and is not calibrated to the benchmark results; the empirical normalized times are compared to the model qualitatively, not used to define the model's parameters. The optimization comparison in Section III B 1 (SLSQP and L-BFGS-B versus Nelder-Mead) is an independent comparison of algorithms using the AD gradients, and no result is assumed in its setup. The paper includes self-citations to the authors' own AWAKE beamline design [2] and distribution-tailoring work [3], but these are used to select a realistic test case and to cite related work, not as load-bearing justification for the efficiency conclusions. The manuscript's own Appendix B explicitly flags a limitation of AD-REV memory scaling to deep lattices, which further indicates that the paper is not presenting a result forced by construction. The skeptic concern about a possibly serial CPU finite-difference baseline is a legitimate benchmarking-robustness question, but it is not a circularity: it does not show that any prediction is equivalent to its input by definition. Overall, no circular step can be exhibited, so the circularity score is minimal; the one-point score reflects only the presence of minor, non-load-bearing self-citations.
Assumptions & free parameters
free parameters (1)
- Central finite-difference step size epsilon =
1e-6
assumptions (4)
- domain assumption Transfer maps are differentiable with respect to parameters; non-differentiable effects (fringe fields, particle losses) are neglected.
- standard math The computational cost model in Eq. (11) with kappa_dp ~ O(1) and constant per-operation cost cma captures the dominant behavior.
- domain assumption JAX reverse-mode AD correctly computes gradients and its implementation overhead is captured by the model.
- domain assumption The benchmark hardware (Intel Xeon Silver 4110, Tesla V100) is representative of resources available to accelerator-physics practitioners.
Cite this review
Pith. "Pith review of Differentiable simulations for particle tracking in accelerators: analysis, benchmarking and optimization." pith.science (2026). https://pith.science/paper/MCV562U2
@misc{pith2026250708476,
author = {Pith},
title = {Pith review of: Differentiable simulations for particle tracking in accelerators: analysis, benchmarking and optimization},
year = {2026},
howpublished = {\url{https://pith.science/paper/MCV562U2}},
note = {Machine review of arXiv:2507.08476}
}
read the original abstract
Optimization of beamlines and lattices is a common problem in accelerator physics, which is usually solved with semi-analytical methods and numerical optimization routines. However, these are usually of the gradient-free or finite-differences type, whose computational cost grows quickly with the number of optimization parameters. On the other hand, the cost of gradient-based optimization can scale well with the number of parameters, but only if the computation of the gradient is itself efficient (e.g. not via finite differences, which are inefficient). Recently, there has been an emergence of so-called "differentiable" codes that efficiently provide gradients. Nevertheless, analysis and benchmarking comparisons of these techniques have largely been absent from the literature. In this work, we develop our own differentiable code, via auto-differentiation. We analyze and benchmark differentiability against finite differences in two test cases, a space-charge FODO cell from the literature and a realistic future beamline at CERN. The analyses and benchmarking are done both theoretically and empirically. Finally, we embed the gradient provided by the differentiable code in gradient-based optimization routines and compare with gradient-free methods. This work offers the first such analysis and benchmarking and thus contributes towards the development of more efficient and performant particle accelerators.
Figures
Figures from the paper (8 more)
Reference graph
Works this paper leans on
-
[1]
It is a 6D tracking prob- lem of 5000 particles over a FODO cell with space-charge effects[26]
Single cell This case is taken from [ 18]. It is a 6D tracking prob- lem of 5000 particles over a FODO cell with space-charge effects[26]. Following the reference, each of the drifts is sliced into segments of 0.05 m and the quadrupoles 4 are sliced into segments of 0.025 m. The drift segments are integrated via drift-kick-drift, with the kick due to spac...
-
[2]
Multiple cells We slightly modify the FODO cell of the previous sec- tion. To slightly reduce the computational cost of each cell, we reduce the number of integration steps of a half slice from 20 to 5. To avoid two consecutive drifts, we dou- ble the length of the first drift and remove the last drift. Taking this modified cell, we make a beamline of up ...
work page 2000
-
[3]
4: Ratio of execution times: FD over AD – GPU
Optimization Even if AD is more efficient than FD at computing the gradient, it does not necessarily follow that gradient-based optimization with AD is more efficient than gradient-free 6 10002000300040005000600070008000900010000 number of particles 1 2 3 4 5 number of cells 0.32 0.5 0.61 0.69 0.69 0.81 1.1 1.2 1.2 1.4 0.33 0.58 0.78 0.81 0.85 1.1 1.6 1.7...
-
[4]
Example As explained in Section II C, auto-differentiation rep- resents the numerous computational operations of a computationally-complex function as directed acyclic graph (DAG) of, usually unary or binary, elementary operations. Figure 11 shows such a graph for the hor- izontal kick of a thin sextupole, which is given by the expression ∆x′ = k2 2 (y2 −...
-
[5]
Cost analysis The AD framework also presents an alternative way to derive the computational cost of differentiation. Here, as in most computational machines, we consider a DAG of unary or binary operations, although the same could be done for any arity. The first step is to compute the primal via a forward pass, which costs CP = nX i=1 c(i) p , (A14) wher...
work page 2000
-
[6]
Iselin, The mad program, in Computing in Accelerator Design and Operation , edited by W
C. Iselin, The mad program, in Computing in Accelerator Design and Operation , edited by W. Busse and R. Zelazny (Springer Berlin Heidelberg, Berlin, Heidelberg, 1984)
work page 1984
-
[7]
R. Ramjiawan, S. D¨ obert, J. Farmer, E. Gschwendtner, F. M. Velotti, L. Verra, G. Z. Della Porta, V. Bencini, and P. N. Burrows, Design and operation of transfer lines for plasma wakefield accelerators using numerical optimizers, Phys. Rev. Accel. Beams 25, 101602 (2022)
work page 2022
-
[8]
M. Remta, F. Velotti, F. Huhn, P. Arrutia Sota, and Y. Dutheil, Towards tailored beam distri- butions for fixed target experiments at cern, in Proc. 16th International Particle Accelerator Conference , IPAC’25 - 16th International Particle Accelerator Confer- ence No. 16 (JACoW Publishing, Geneva, Switzerland,
Show all 40 references
-
[9]
L. Yang, D. Robin, F. Sannibale, C. Steier, and W. Wan, Global optimization of an accelerator lattice using mul- tiobjective genetic algorithms, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 609, ...
2009
-
[10]
Deniau et al., MethodicalAcceleratorDesign/MAD-X: 5.09.01 (2023)
L. Deniau et al., MethodicalAcceleratorDesign/MAD-X: 5.09.01 (2023)
2023
-
[11]
J. A. Nelder and R. Mead, A simplex method for function minimization, The Computer Journal 7, 308 (1965)
1965
-
[12]
Iadarola et al., Xsuite: An Integrated Beam Physics Simulation Framework, JACoW HB2023, TUA2I1 (2024), arXiv:2310.00317 [physics.acc-ph]
G. Iadarola et al., Xsuite: An Integrated Beam Physics Simulation Framework, JACoW HB2023, TUA2I1 (2024), arXiv:2310.00317 [physics.acc-ph]
2024 arXiv
-
[13]
Hansen, A
N. Hansen, A. Auger, R. Ros, S. Finck, and P. Poˇ s ´ ık, Comparing results of 31 algorithms from the black-box optimization benchmarking bbob-2009 (Association for Computing Machinery, New York, NY, USA, 2010) p. 1689–1696
2009
-
[14]
What is considered mid and high dimensional depends on the problem and the computational resources available. 10
-
[15]
Jameson, Aerodynamic design via control theory, Jour- nal of Scientific Computing 3, 233 (1988)
A. Jameson, Aerodynamic design via control theory, Jour- nal of Scientific Computing 3, 233 (1988)
1988
-
[16]
Magri, Adjoint methods as design tools in thermoa- coustics, Applied Mechanics Reviews 71, 020801 (2019)
L. Magri, Adjoint methods as design tools in thermoa- coustics, Applied Mechanics Reviews 71, 020801 (2019)
2019
-
[17]
Our World in Data, Computation used to train notable ai systems, by domain (2024)
2024
-
[18]
Ansel, E
J. Ansel, E. Yang, H. He, N. Gimelshein, A. Jain, M. Voznesensky, B. Bao, P. Bell, D. Berard, E. Burovski, G. Chauhan, A. Chourdia, W. Constable, A. Desmaison, Z. DeVito, E. Ellison, W. Feng, J. Gong, M. Gschwind, B. Hirsh, S. Huang, K. Kalambarkar, L. Kirsch, M. La- zos, M. L...
2024
-
[19]
Abadi, A
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, R. Jozefowicz, Y. Jia, L. Kaiser, M. Kudlur, J. Levenberg, D. Man´ e, M. Schuster, R. Monga, S. Moore, D. Mur- r...
2015
-
[20]
Bradbury, R
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. Van- derPlas, S. Wanderman-Milne, and Q. Zhang, JAX: com- posable transformations of Python+NumPy programs (2025)
2025
-
[21]
Gonzalez-Aguilera, Y.-K
J. Gonzalez-Aguilera, Y.-K. Kim, R. Rous- sel, A. Edelen, and C. Mayes, Towards fully differentiable accelerator modeling, in Proc. 14th International Particle Accelerator Conference , IPAC’23 - 14th International Particle Accelerator Confer- ence No. 14 (JACoW Publishing, Gen...
-
[22]
Kaiser, C
J. Kaiser, C. Xu, A. Eichler, and A. Santamaria Garcia, Bridging the gap between machine learning and particle accelerator physics with high-speed, differentiable simula- tions, Phys. Rev. Accel. Beams 27, 054601 (2024)
2024
-
[23]
J. Wan, H. Alamprese, C. Ratcliff, J. Qiang, and Y. Hao, JuTrack: a Julia package for auto-differentiable accelerator modeling and particle tracking (2024), arXiv:2409.20522 [physics.acc-ph]
2024 arXiv
-
[24]
Deniau, Mad-ng, a standalone multiplatform tool for linear and non-linear optics design and optimisation (2025), arXiv:2412.16006 [cs.CE]
L. Deniau, Mad-ng, a standalone multiplatform tool for linear and non-linear optics design and optimisation (2025), arXiv:2412.16006 [cs.CE]
2025 arXiv
-
[25]
J. Wan, H. Alamprese, C. Ratcliff, J. Qiang, and Y. Hao, JuTrack code repository (2025), Accessed: 2025-06-01
2025
-
[26]
Not to be confused with forward FD, which is unrelated
-
[27]
CD = 2CP
For example, if the primal is f (x, y) = xy, then ∂pf = y ∂px + x ∂py is y and x, i.e. CD = 2CP
-
[28]
Practical limitations, like memory, aside
-
[29]
A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of Machine Learning Research 18, 1 (2018)
2018
-
[30]
The code used in this paper can be found here
-
[31]
Details can be found in the reference
-
[32]
The reason for the slight discrepancies is unknown
-
[33]
This estimate is obtained via the func.lower(x).compile().cost analysis() method in JAX
-
[34]
Coherent synchrotron radiation is the only effect that could be considered but for the purpose of this paper, it was considered negligible
-
[35]
The order of the elements and the values of the param- eters of the initial design can be found in the project’s repository
-
[36]
2a vs Fig
Confirmed by the primals’ execution times (see Fig. 2a vs Fig. 8a)
-
[37]
Kraft, Algorithm 733: Tomp–fortran modules for op- timal control calculations, ACM Transactions on Mathe- matical Software 20, 262–281 (1994)
D. Kraft, Algorithm 733: Tomp–fortran modules for op- timal control calculations, ACM Transactions on Mathe- matical Software 20, 262–281 (1994)
1994
-
[38]
C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound- constrained optimization, ACM Trans. Math. Softw. 23, 550–560 (1997)
1997
-
[39]
It is perhaps then not surprising that one of the few publications [ 35] on application of AD to beamlines is on optimizing a neural network that generates the ini- tial distribution, and not on optimizing the beam-line elements
-
[40]
Roussel, J
R. Roussel, J. P. Gonzalez-Aguilera, E. Wisniewski, A. Ody, W. Liu, J. Power, Y.-K. Kim, and A. Edelen, Effi- cient six-dimensional phase space reconstructions from ex- perimental measurements using generative machine learn- ing, Phys. Rev. Accel. Beams 27, 094601 (2024)
2024
Reviewed August 6, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.