Pith. sign in

REVIEW 3 major objections 5 minor 24 references

A parallel algorithm for generating Pareto-optimal radiosurgery treatment plans

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

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

Pith's one-line read A GPU-parallel ADMM solver can produce hundreds of Pareto-optimal Gamma Knife plans in minutes, matching clinical simplex quality.

desk verdict Solid engineering contribution: ADMM parallelized over weight vectors gives clinically matching plan quality at large speedups, but the 'Pareto-optimal' label is loose and the empirical tuning rules deserve a closer look. read the letter →

arxiv 2509.08602 v1 pith:SQFWLPKB submitted 2025-09-10 physics.med-ph math.OC

Joakim da Silva , Daniel Hernández Escobar , Tor Kjellsson Lindblom , H{aa}kan Nordström , Jens Sjölund This is my paper · ORCID
classification physics.med-phmath.OC
keywords GammaKniferadiosurgeryinversetreatmentplanningmulticriteriaoptimizationParetonavigationADMMGPUparallelcomputinglinearprogrammingtwo-passdose
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

Radiosurgery planning today is iterative: the planner changes a slider, reruns the optimizer, and checks the plan. The paper aims to remove that waiting by solving many trade-off variants at once with a GPU-parallel version of ADMM, the alternating direction method of multipliers, applied to the linear program already used in clinical Gamma Knife planning. The key structural observation is that clinical trade-off weights enter only through the bound constraints, so after one cached matrix factorization each iteration becomes dense matrix products plus elementwise projections, and hundreds of weight vectors can be advanced together. On 20 clinical cases the authors report batch generation of up to 441 plans in 3–110 seconds on a GPU—under two minutes even for the largest targets—with clinical metrics matching the single-threaded simplex solver on 3600 comparison runs. If it holds, this enables interactive Pareto navigation of clinical trade-offs on the day of treatment.

What carries the argument

The workhorse is a variable-splitting ADMM in which the original LP's bound constraints are projected elementwise while the equality-constrained step reduces to a cached Schur complement, the matrix AA^T+I, whose factorization turns each iteration into dense matrix products plus clipping; stacking weight vectors as matrix columns makes a whole ensemble of plans one batched computation. A two-pass scheme with an overlapping low-dose point sampler couples the second-pass objective to the first-pass solution without breaking parallelism.

What would settle it

Take a previously unseen case with larger total target volume or more targets than the tested set, run the same fixed 3000-iteration ADMM over the full 21x21 slider grid, and compare coverage, selectivity, gradient index, and beam-on time against the simplex solver; a systematic gap larger than the reported few percent in any metric would show the tuned regime does not generalize.

Watch

Extended reading notes

Core claim

The central claim is that a GPU-parallel ADMM solver, specialized to the dual linear-program formulation of Gamma Knife inverse planning, can generate several hundred Pareto-optimal treatment plans (up to 441, corresponding to a dense grid of clinical slider settings) in 3–110 seconds across the tested cases, and in under two minutes even for the largest targets. Compared with sequentially running the single-threaded dual-simplex solver, the authors report speedups of 54–1500 times on GPU and 1.6–97 on CPU for batches of 81 plans, with the largest speedups on the hardest cases. Plan quality was checked by rerunning both solvers 20 times per case and weight combination on 20 cases (3600 solve

Load-bearing premise

The fixed tuning rules for the ADMM step size and column scaling, and the fixed 3000-iteration budget, were chosen on 10 validation cases and are assumed to hold for any new anatomy and any slider combination; a case outside that tuned regime could produce plans whose quality deviates from the clinical reference.

Editorial extensions

If this is right

  • A full 21x21 grid of slider settings, 441 plans, fits within the time frame of a clinical session on one GPU, even for large targets.
  • Plan quality and run-to-run variability match the simplex-based clinical reference within a few percent across the four standard clinical metrics.
  • Interactive, interruption-free Pareto navigation becomes practical: clinicians can explore trade-offs in real time rather than rerun the optimizer per slider move.
  • The CPU implementation also gives 1.6–97x speedups over single-threaded simplex for 81 plans, so centers without GPUs still benefit.

Reading between the lines

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

  • The authors tuned the step size and column scaling on 10 validation cases and fixed them; an obvious extension is to make these parameters adaptive per weight vector or derived from the constraint matrix so the method does not rely on a pre-tuned regime.
  • The paper compares against sequential simplex; in a realistic clinical setting where several plans are wanted, memory limits prevented more than four parallel simplex solves for large cases, so the practical speed advantage of batch ADMM is likely larger than the headline simplex ratio suggests.
  • The same batch-parallel pattern should transfer to other convex fluence-map optimization problems where objectives are piecewise-linearized; the paper states this expectation without demonstrating it, so testing it on IMRT-style problems would be a natural next step.
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.

Referee Report

3 major / 5 minor

Summary. The paper proposes a parallel ADMM-based solver for the linear-program formulation used in Gamma Knife radiosurgery planning (the Lightning module). The authors derive an ADMM splitting of the dual LP, cache the Schur complement, and parallelize the iterations across many objective-weight vectors to generate a large batch of plans simultaneously. A two-pass low-dose point sampling scheme is introduced to handle gradient-index control in parallel. The method is implemented in Matlab on CPU and GPU and evaluated on 20 clinical test cases (three indications, target volumes 0.66–52 cm³) against Matlab's single-threaded dual-simplex solver. The paper reports speedups of 1.6–97× (CPU) and 54–1500× (GPU) for 81 weight vectors, with 441 plans generated in 3–110 s on GPU, and clinical metrics (coverage, selectivity, gradient index, beam-on time) statistically matching those of the simplex reference across 20 random-seed reruns. The central claim is that the method enables interactive Pareto navigation by generating hundreds of clinically acceptable plans within minutes.

Significance. If the claims hold, this is a practically important contribution: it demonstrates a concrete route to multicriteria navigation for Gamma Knife radiosurgery without per-plan optimization. The evaluation is unusually thorough for a methods paper: a separate validation set is used to tune hyperparameters, 20 test cases span a wide range of target complexity, and clinical metrics are compared over 20 random-seed reruns. The mathematical derivation of the ADMM updates and the Schur-complement caching is sound, and the idea of parallelizing across weight vectors is well matched to the clinical slider-based workflow. However, the paper's central 'Pareto-optimal' and 'without compromising quality' claims rely on empirically tuned hyperparameters and a fixed 3000-iteration budget, and the generalization of these choices beyond the studied cohort is not established. The lack of released code/data further limits independent verification, though the algorithmic description is sufficiently clear to be reimplemented.

major comments (3)
  1. [Section III.C and Section III.D (Eqs. 9–10, fixed 3000 iterations)] The method's robustness rests on the empirical tuning rules for ρ and β (Eqs. 9 and 10) and on the fixed 3000-iteration budget, all chosen using a 10-case validation set and then fixed for the 20 test cases. This is a legitimate hold-out evaluation, but it does not establish that these choices will remain adequate for anatomies outside the studied range (e.g., very different dose-point/isocenter ratios, unusual dose-rate kernels, or different OAR configurations). The authors themselves note that the BOT term was 'somewhat worse' for the largest validation case until β scaling was added, demonstrating sensitivity to problem size. Because the paper presents the method as a turnkey clinical replacement, I ask for: (i) a residual- or duality-gap-based stopping criterion, or at least a per-case convergence check; (ii) objective-suboptimality reporting on the 20 test cases (the 1–3% objective-
  2. [Title, Section I, Section II.C] The title and abstract claim 'Pareto-optimal treatment plans', but the paper does not establish exact Pareto optimality. Exact solutions of the weighted-sum LP with positive weights are Pareto-optimal for the sampled problem, but ADMM with a fixed 3000-iteration budget only gives approximate solutions, and the two-pass low-dose sampling (Section II.C) changes the optimization problem between passes by adapting the low-dose points to the first-pass dose distribution. The final plans are therefore at best approximately Pareto-optimal for an approximation of the clinical objectives. Please qualify the claim (e.g., 'approximately Pareto-optimal' or 'numerically Pareto-optimal'), define the precise Pareto-optimality notion, and report a direct optimality-gap measure for the test cases rather than only clinical metrics, which can be insensitive to small objective deviations.
  3. [Section III.E, Figures 4 and 5] The headline speedups compare a batch-parallel approximate solver (ADMM across many weight vectors on GPU/CPU) to a sequential, single-threaded exact simplex solver. This is disclosed in the text, and the memory argument against parallelizing simplex is reasonable. However, the abstract and conclusions present the speedup as an algorithmic property when it is substantially a batch-throughput and parallelism gain. I recommend reporting also single-instance ADMM times, strong/weak scaling of the batch, and, if possible, a comparison with a multi-threaded or warm-started simplex baseline, or explicitly reframing the numbers as 'batch throughput speedup for interactive navigation'. This would make the central practical claim more precise.
minor comments (5)
  1. [Section III.C, Eqs. (9)–(10)] The symbols N_TS, N_TI, N_R, N_1,LD, N_2,LD are used in Eq. (10) but defined only later in Appendix A. Please add a forward reference or brief definition at first use.
  2. [Section II.C.2, Algorithm 1] The notation in Algorithm 1 is confusing: p_tilde_j^union is 'randomly select ω_j |V_{j-1}^union| points from p_{j-1}^union', but p_{j-1}^union is a set of points, not a volume, so the intended sampling mechanism and the relationship between |V| and |p| should be stated explicitly. Also clarify whether the selection is with or without replacement.
  3. [Section III.A] Typo: 'Intel Xenon Gold' should be 'Intel Xeon Gold'.
  4. [Section III.A and III.C] The GPU implementation uses an explicit inverse of the Schur complement in single precision to avoid repeated decomposition. The paper says this 'could worsen numerical stability' but was not observed as a problem. Please report the condition number of S (or a related diagnostic) for the largest test cases to support this claim.
  5. [Reproducibility] No code or data availability statement is provided. Given that the tuning rules and random sampling are central to the method, a statement on availability (even for the anonymized test cases) would strengthen reproducibility.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the ADMM method is benchmarked against an external simplex solver on held-out test cases; hyperparameters are tuned on a separate validation set.

full rationale

The paper's central claim is that its GPU/CPU ADMM implementation produces plans of comparable quality to the simplex-based clinical solver in far less time. This is an empirical comparison, not a derivation from a fitted quantity. The hyperparameters ρ (Eq. 9) and β (Eq. 10) were tuned on a separate 10-case validation set and then fixed for the 20-case test set; none of the reported clinical metrics or timings are used to fit those parameters. The LP formulation is inherited from the authors' prior work (Sjölund et al. 2019), but that is an input modeling choice, not a predicted outcome; the paper's contribution is the parallel ADMM adaptation and overlapping low-dose sampling, both evaluated against an external solver rather than against the paper's own assumptions. The two-pass optimization depends on first-pass solutions only as an iterative sampling heuristic, not as a self-referential prediction. The timing caveat that simplex is single-threaded is a fairness limitation, not a circularity. No equation in the paper defines the quantity it claims to predict, and no success metric reduces to a fitted parameter by construction. Therefore no circularity is present.

Assumptions & free parameters 4 free parameters · 6 assumptions · 0 invented entities

The method introduces no new physical entities. It relies on standard convex optimization assumptions and the existing clinical LP formulation. The main empirical elements are the tuned hyperparameters (rho, beta, iteration count) and the two-pass heuristic, which are all calibrated on validation data rather than derived from first principles.

free parameters (4)
  • ADMM step size coefficient (rho rule) = 2.5e-3
    The coefficient in Equation (9) was selected via grid search on a validation set to match simplex cost terms. It is a free parameter tuned to data.
  • Column scaling denominator (beta rule) = 2000
    The denominator in Equation (10) was chosen based on validation cases to scale BOT columns. It is a hand-tuned constant.
  • Number of ADMM iterations = 3000 per pass
    Chosen as a fixed iteration count after observing convergence on validation cases; not derived from a stopping criterion.
  • Low-dose threshold interval Delta = Not specified in text
    The interval [D_LD, D_LD + Delta) defining low-dose volumes is referenced but its value is not given; presumably inherited from the clinical software.
assumptions (6)
  • domain assumption The dual LP (Equation 1) is equivalent to the primal multicriteria problem of Sjölund et al.
    The paper relies on the duality formulation from reference [1] without re-deriving it. This is a standard result but the specific mapping of weights to bounds is taken from prior work.
  • standard math ADMM converges to the solution of the LP for the chosen splitting and updates.
    ADMM convergence is assumed from Boyd et al. [18]; the paper does not provide a bespoke convergence proof for this problem.
  • domain assumption Isocenter placement is independent of objective weights.
    The paper uses the Lightning isocenter set, which is a function only of target geometry. This is stated in Section II.A and is load-bearing for solving all weight vectors in parallel.
  • domain assumption Sampling dose points from target surfaces, rings, and low-dose regions is sufficient to capture plan quality.
    The paper follows the sampling strategy of Sjölund et al. and states that more sophisticated strategies would not change conclusions. This is a modeling assumption.
  • domain assumption The two-pass optimization heuristic yields clinically acceptable plans for all weight vectors.
    The adaptive low-dose volume selection is a heuristic borrowed from clinical practice; the paper does not prove optimality, only empirical equivalence to simplex on test cases.
  • domain assumption OAR max-dose constraints can be enforced by sampling points on OAR surfaces.
    This relies on the assumption that no isocenter lies inside an OAR, stated in Section II.A. It is a geometric assumption that may fail in some anatomies.

how reviews work

0 comments
Cite this review

Pith. "Pith review of A parallel algorithm for generating Pareto-optimal radiosurgery treatment plans." pith.science (2026). https://pith.science/paper/SQFWLPKB

@misc{pith2026250908602,
  author       = {Pith},
  title        = {Pith review of: A parallel algorithm for generating Pareto-optimal radiosurgery treatment plans},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/SQFWLPKB}},
  note         = {Machine review of arXiv:2509.08602}
}
read the original abstract

Using inverse planning tools to create radiotherapy treatment plans is an iterative process, where clinical trade-offs are explored by changing the relative importance of different objectives and rerunning the optimizer until a desirable plan is found. We seek to optimize hundreds of radiosurgery treatment plans, corresponding to different weightings of objectives, fast enough to incorporate interactive Pareto navigation of clinical trade-offs into the clinical workflow. We apply the alternating direction method of multipliers (ADMM) to the linear-program formulation of the optimization problem used in Lightning. We implement both a CPU and a GPU version of ADMM in Matlab and compare them to Matlab's built-in, single-threaded dual-simplex solver. The ADMM implementation is adapted to the optimization procedure used in the clinical software, with a bespoke algorithm for maximizing overlap between low-dose points for different objective weights. The method is evaluated on a test dataset consisting of 20 cases from three different indications, with between one and nine targets and total target volumes ranging from 0.66 to 52 cm3, yielding speedups of 1.6-97 and 54-1500 times on CPU and GPU, respectively, compared to simplex. Plan quality was evaluated by rerunning the ADMM optimization 20 times, each with a different random seed, for each test case and for nine objective weightings per case. The resulting clinical metrics closely mimicked those obtained when rerunning the simplex solver, verifying the validity of the method. In conclusion, we show how ADMM can be adapted for radiosurgery plan optimization, allowing hundreds of high-quality Gamma Knife treatment plans to be created in under two minutes on a single GPU, also for very large cases.

Figures

Figures reproduced from arXiv: 2509.08602 by the authors.

Figure 1
Figure 1. Distributions of clinical metrics for case ME10 after the second optimization pass. [PITH_FULL_IMAGE:figures/full_fig_p012_1.png] view at source ↗
Figure 2
Figure 2. Distributions of clinical metrics after the second optimization pass, normalized to [PITH_FULL_IMAGE:figures/full_fig_p013_2.png] view at source ↗
Figure 3
Figure 3. Distributions of clinical metrics after the second optimization pass, normalized [PITH_FULL_IMAGE:figures/full_fig_p014_3.png] view at source ↗
Figures from the paper (2 more)
Figure 4
Figure 4. Figure 4: Log-log plots of ADMM optimization times in seconds for two passes of 3000 [PITH_FULL_IMAGE:figures/full_fig_p016_4.png]
Figure 5
Figure 5. Figure 5: Log-log plots of speedups when using the ADMM method on GPU (rings) and CPU [PITH_FULL_IMAGE:figures/full_fig_p017_5.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

24 extracted references · 1 linked inside Pith

  1. [1]

    o lund, S. Riad, M. Hennix, and H. Nordstr\

    J. Sj\" o lund, S. Riad, M. Hennix, and H. Nordstr\" o m, A linear programming approach to inverse planning in Gamma Knife radiosurgery, Medical physics 46 , 1533--1544 (2019)

    2019

  2. [2]

    Ham \'a c kov \'a , J

    L. Ham \'a c kov \'a , J. Novotný Jr, M. Farn \' kov \'a , R. Li s c \'a k, G. S imonov \'a , D. Urgo s \' k, and M. Schmitt, Comparative study evaluating a new dose optimization software for Gamma Knife treatment planning--comparison of 80 challenging treatment plans, Stereotactic and Functional Neurosurgery (2025)

    2025

  3. [3]

    D. L. Craft, T. F. Halabi, H. A. Shih, and T. R. Bortfeld, Approximating convex Pareto surfaces in multiobjective radiotherapy planning, Medical physics 33 , 3399--3407 (2006)

    2006

  4. [4]

    H. E. Romeijn, J. F. Dempsey, and J. G. Li, A unifying framework for multi-criteria fluence map optimization models, Physics in Medicine & Biology 49 , 1991 (2004)

    1991

  5. [5]

    Breedveld, P

    S. Breedveld, P. R. Storchi, P. W. Voet, and B. J. Heijmen, iCycle : Integrated, multicriterial beam angle, and profile optimization for generation of coplanar and noncoplanar IMRT plans, Medical physics 39 , 951--963 (2012)

    2012

  6. [6]

    Tolakanahalli, D

    R. Tolakanahalli, D. J. J. Wieczorek, Y. C. Lee, M. C. Tom, M. D. Hall, M. W. McDermott, M. P. Mehta, R. Kotecha, and A. N. Gutierrez, OptImal Gamma kNife lIghTnIng sOlutioN (IGNITION) score to characterize the solution space of the Gamma Knife FIP optimizer for stereotactic radiosurgery, Journal of applied clinical medical physics 24 , e13936 (2023)

    2023

  7. [7]

    Miettinen, F

    K. Miettinen, F. Ruiz, and A. P. Wierzbicki, Introduction to multiobjective optimization: interactive approaches, in Multiobjective optimization: interactive and evolutionary approaches , pages 27--57, Springer, 2008

    2008

  8. [8]

    Monz, K.-H

    M. Monz, K.-H. K \"u fer, T. R. Bortfeld, and C. Thieke, Pareto navigation—algorithmic foundation of interactive multi-criteria IMRT planning, Physics in Medicine & Biology 53 , 985 (2008)

    2008

Show all 24 references
  1. [9]

    Bokrantz and A

    R. Bokrantz and A. Forsgren, An algorithm for approximating convex Pareto surfaces based on dual techniques, INFORMS Journal on Computing 25 , 377--393 (2013)

    2013

  2. [10]

    Glowinski and A

    R. Glowinski and A. Marroco, Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de Dirichlet non linéaires, Revue française d'automatique, informatique, recherche opérationnelle. Analyse numérique 9 , 41--76 (1975)

    1975

  3. [11]

    Gabay and B

    D. Gabay and B. Mercier, A dual algorithm for the solution of nonlinear variational problems via finite element approximation, Computers & mathematics with applications (1987) 2 , 17--40 (1976)

    1987

  4. [12]

    Q. J. Wu, V. Chankong, S. Jitprapaikulsarn, B. W. Wessels, D. B. Einstein, B. Mathayomchan, and T. J. Kinsella, Real-time inverse planning for Gamma Knife radiosurgery, Medical physics 30 , 2988--2995 (2003)

    2003

  5. [13]

    Doudareva, K

    E. Doudareva, K. Ghobadi, D. M. Aleman, M. Ruschin, and D. A. Jaffray, Skeletonization for isocentre selection in Gamma Knife Perfexion , Top 23 , 369--385 (2015)

    2015

  6. [14]

    Ghobadi, H

    K. Ghobadi, H. R. Ghaffari, D. M. Aleman, D. A. Jaffray, and M. Ruschin, Automated treatment planning for a dedicated multi-source intracranial radiosurgery treatment unit using projected gradient and grassfire algorithms, Medical physics 39 , 3134--3141 (2012)

    2012

  7. [15]

    Fountain, K

    L. Fountain, K. Khedriliraviasl, S. Mahmoudzadeh, and H. Mahmoudzadeh, Dose-based constraint generation for large-scale IMRT optimization, INFOR: Information Systems and Operational Research 60 , 1--19 (2022)

    2022

  8. [16]

    S. Mair, A. Fu, and J. Sj \"o lund, Efficient radiation treatment planning based on voxel importance, Physics in Medicine & Biology 69 , 165031 (2024)

    2024

  9. [17]

    Quarz, L

    A. Quarz, L. Volz, C. H. Antink, M. Durante, and C. Graeff, Deep learning-based voxel sampling for particle therapy treatment planning, Physics in Medicine & Biology 69 , 155014 (2024)

    2024

  10. [18]

    Boyd et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine learning 3 , 1--122 (2011)

    S. Boyd et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine learning 3 , 1--122 (2011)

    2011

  11. [19]

    Nair et al., Solving mixed integer programs using neural networks, arXiv preprint arXiv:2012.13349 (2020)

    V. Nair et al., Solving mixed integer programs using neural networks, arXiv preprint arXiv:2012.13349 (2020)

    2012 arXiv

  12. [20]

    Zarepisheh, M

    M. Zarepisheh, M. Shakourifar, G. Trigila, P. Ghomi, S. Couzens, A. Abebe, L. Nore \ n a, W. Shang, S. B. Jiang, and Y. Zinchenko, A moment-based approach for DVH -guided radiotherapy treatment plan optimization, Physics in Medicine & Biology 58 , 1869 (2013)

    2013

  13. [21]

    A. Fu, B. Ungun, L. Xing, and S. Boyd, A convex optimization approach to radiation treatment planning with dose constraints, Optimization and engineering 20 , 277--300 (2019)

    2019

  14. [22]

    , Optimization Toolbox (R2024b) , The MathWorks Inc., Natick, Massachusetts, United States, 2024

    The MathWorks Inc. , Optimization Toolbox (R2024b) , The MathWorks Inc., Natick, Massachusetts, United States, 2024

    2024

  15. [23]

    H. E. Romeijn, R. K. Ahuja, J. F. Dempsey, A. Kumar, and J. G. Li, A novel linear programming approach to fluence map optimization for intensity modulated radiation therapy treatment planning, Physics in Medicine & Biology 48 , 3521 (2003)

    2003

  16. [24]

    D. M. Shepard, M. A. Earl, X. A. Li, S. Naqvi, and C. Yu, Direct aperture optimization: a turnkey solution for step-and-shoot IMRT , Medical physics 29 , 1007--1018 (2002)

    2002

Pith tools

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