From Sequential to Parallel: Reformulating Dynamic Programming as GPU Kernels for Large-Scale Stochastic Combinatorial Optimization
Authors: Jingyi Zhao
Linxin Yang
Haohua ZHANG
Qile He
Tian Ding
Venue: ICLR 2026 Poster (2026)
Keywords: CUDA
GPU computing
Large-scale optimization
Dynamic Programming
Stochastic Optimization
Primary area: optimization

Abstract: Dynamic programming (DP) is central to combinatorial optimization, optimal control, and reinforcement learning, yet its perceived sequentiality has long hindered scalability. We introduce a general-purpose GPU framework that reformulates broad classes of forward DP recursions as batched min--plus matrix--vector products over layered DAGs, collapsing actions into masked state-to-state transitions that map directly to GPU kernels. This approach removes a major bottleneck in scenario-based stochastic programming (SP), where the use of DP has traditionally restricted the number of scenarios due to excessive computational cost.
Our framework exposes massive parallelism across scenarios, transition layers, and, when applicable, route or action options, via self-designed GPU kernels that implement Bellman updates with warp-/block-level reductions and numerically safe masking. In a single GPU pass, these kernels can process over $10^6$ uncertainty realizations, far beyond the capacity of prior scenario-based methods. We demonstrate the approach in two canonical SP applications: (i) a vectorized split operator for the capacitated vehicle routing problem with stochastic demand, exploiting **2D** parallelism (scenarios $\times$ transitions); and (ii) a forward inventory reinsertion DP under an order-up-to policy, exploiting **3D** parallelism (scenarios $\times$ inventory transitions $\times$ route options). Across benchmarks, the implementation scales nearly linearly in the number of scenarios and achieves one to three orders of magnitude speedups over multithreaded CPU baselines, yielding tighter SAA estimates and consistently stronger first-stage decisions under identical wall-clock budgets. Viewed as hardware-aware software primitives, our min--plus DP kernels offer a drop-in path to scalable, GPU-accelerated stochastic discrete optimization.

Published as a conference paper at ICLR 2026
FROM SEQUENTIAL TO PARALLEL: REFORMULAT- ING DYNAMIC PROGRAMMING AS GPU KERNELS FOR LARGE-SCALE STOCHASTIC COMBINATORIAL OPTI-
MIZATION
Jingyi Zhao', Linxin Yang!2, Haohua Zhang', Qile Hc², Tian Ding Shenzhen Research Institute Big Data, Shenzhen, China
2. School of Data Science, The Chinese University of Hongkong, Shenzhen, China
Shenzhen International Center for Industrial and Applied Mathematics, Shenzhen, China
4. AutoKemne,, Shenzhen, China
ABSTRACT
A major bottleneck in scenario based Sample Average Approximation (SAA) for stochastic programming (SP) is the cost of solving an exact second-stage problem for every scenario, especially when each scenario contains NP-hard combina- torial structure. This has led much of the SP literature to restrict the second stage linear or simplified models. We develop GPU-based framework that makes full-fidelity integer second-stage models tractable at scale. The key innovation is a set of hardware-aware scenario batched GPU kernels that expose parallelism across scenarios, dynamic-programmit programming (DP) layers, and route or action options, enabling Bellman updates be executed single pass over more than 10 real- izations. We evaluate the approach in two representative SP settings: a vectorized split operator for stochastic vehicle routing and a DP for inventory reinsertion. Implementation scales nearly linearly in the number of scenarios and achieves one-two to four-five orders of magnitude speedup, allowing far larger scenario sets and reliably stronger first-stage decisions. The computational leverage di- rectly improves decision quality: much larger scenario sets and many more first- stage candidates can be evaluated within fixed time budgets, consistently yielding stronger SAA solutions. Our results show that full-fidelity integer second-stage models are tractable at scales previously considered impossible, providing prac- tical path to large-scale, realistic stochastic discrete optimization.
INTRODUCTION
A central challenge in scenario-based stochastie programming (SP) is the computational burden of solving an exact second-stage problem for every scenario, especially when each scenario embeds an NP-hard combinatorial structure. This difficulty has led much of the SP literature to restrict the second stage to linear programs or simplified models (Birge & Louveaux] 1997] Shapiro et al., 202T,, sacrificing realism and often degrading first-stage decision quality. Our goal to make full- fidelity, integer second-stage models computationally viable at scale, without relying on surrogate learning models or structural relaxations.
A key observation is that once a first-stage decision is fixed, second recourse evaluations across scenarios are independent. Many of these evaluations rely on dynamic programming (DP), a foun- dational paradigm in optimization and control (Bertsekas) (012) and a core ingredient in both exact and heuristic combinatorial optimization. DP underlies classical algorithms such as the Held- Karp procedure for TSP (Held & Karp 1971). pseudo-polynomic knapsack solvers (Shapiro 1968), and
shortest path algorithms (Bellman 1958 Dijkstra 1959). It appears as a subroutine in vehicle and inventory routing via split operators (Prins 20044 Vidal 2016) and resource- con.urcceconttriieeeddeeeeeeeedd shortest- path pricing (Feillet et al.) 2004 Imnich & Desaulniers 2005) and large-scale integer programs (Barhhart et al. 1998). DP also plays central role in metaheuristics (Vidal et al.] 2012 Zhao al.

---

Published as a conference paper at ICLR 2026
2022
and multistage stochastic optimization (Powell 2007 Bertsekas & Tsitsiklis 1995 Staddon
2020
Despite this ubiquity, most DP implementations are inherently sequential, obscuring the substantial parallelism available across scenarios and within DP transitions. While recent progress GPU acceleration has benefited continuous optimization (Schubiger 2019: Lu & Yang 2023 Bishop
ett al. 2024 Liu et al. 2024), no comparable framework exists for discrete, combinatorial, scenario- based problems.
We address this gap by developing a GPU-based framework that executes second-stage dynamic programs in scemario-batched and multidimensional manner. The key novelty lies in custom GPU kernels that expose parallelism simultaneously across scenarios, DP layers, and route or action op- tions. These kernels implement Bellman updates using warp-/block-leve reductions with numeri- cally safe masking, enabling single GPU pass to evaluate over 10t uncertainty realizations far beyond what existing SP approaches can handle.
We demonstrate the approach on two benchmark problems that capture core challenges in stochastic routing and inventory optimization: (i) a vectorized split operator for capacitated vehicle routing with stochastic demand (CVRPSD), and (ii) an inventory reinsertion DP un- der an order-up-to policy. Across benchmarks, our implementation scales nearly linearly with the number of scenarios and achieves substantial speedups over CPU baselines- one-two orders of magnitude for CVRPSD to four-five orders for DSIRP. Additional experi- ments show that these gains translate directly into better decision quality: the ability to eval- uate far more scenarios and far more first-stage candidates within the same time budget con- sistently yields stronger solutions. Our implementation is available htttps://githh: and https://
github wtt.
In summary, our main contributions are as follows.
1. We expose and formalize the multidimensional parallelism inherent in second-stage dy- namic programs for stochastic combinatorial optimization, and demonstrate that this struc- ture enables scenario-batcher exccution at scales previously out of reach, surpassing 10t scenarios.
2. We design hardware-aware GPU kernels that exploit concurrency across scenarios, DP layers, and route/action options. These kernels implement efficient Bellman reductions with numerically safe masking, enabling high-throughput evaluation of massive scenario batches.
3. We present extensive empirical evidence, including comparisons against the extensive-for form MILP solved by Gurobi for CVRPSD and against the state-of-the-art DSIRP algorithm (Coelho al., 2012). Our method remains tractable for far larger scenario and fre- quently obtains higher-quality solutions.
4. We demonstrate that GPU-accelerated second-stage evaluation substantially strengthens metaheuristic search: orders-of -magnnitue more first-stage candidates and far larger see- nario sets can be explored within the same time budget. SAA, solution quality improves reliably as the number of scenarios increases, making computational breakthrough di- rectly valuable for obtaining high-quality stochastic solutions.
5. We provide a detailed GPU-memory study showing that realistic problem instances with up 10 scenarios fit within standard 1IGB GPU, confirming that memory is not the limiting factor at scale. We also outline a general recipe for converting DP subroutines into tigh-throughput GPU primitives.
2 GENERIC DYNAMIC PROGRAMMING FRAMEWORK
Preliminary Knowledge. We consider finite-horizon dynamic program over stages 1,...,T, T, starting from an initial state $1 E Sy. At each stage t. the system is in state s, E Sü, an action a, Ar(st) is chosen, and the system moves 8t+1 either deterministically via go(st, or stochastically according to (5 St, w). Each transition incurs stage cost (s,, under scenario
2

---

Published as a conference paper at ICLR 2026
Let (s) denote the minimum cumulative cost to reach St at stage t. The recursion initializes as J ($1) and evolves by
J (++1(0)=
min sES. (EA)(s) si(s,a)ms'
(JF(s) + (s, Aal E S4+1.
(1)
The objective is to minimize the expected terminal cost E min [s)]
2.1 TRANSITION-BASER FORMULATION.
To enable GPU-friendLy computation, express the recursion in terms of state-to-state transition costs. For each stage and scenario w, define
AP (s, (s),
inf (s,a), a (EAr(a) ge ==,,a)==
(s, +oo no feasible transition exists.
The forward recursion then becomes
18'1 (sss sES,
This formulation results in stage-wise min-plus update that can be evaluated independently for each stage, making well aligned with batched GPU execution.
2.2 MIN-PLUS MATRIX FORMULATION.
Let S, = (1,...,m,, m,) index the state space. Define:
AT (i,j) A i, j)= J e Rm.
Then, the Bellman update becomes a matrix-vector product in the (min, +) semiring:
J = (A7) 2 J = (nin(4 + JF
(2)
This min-plus formulation enables efficient GPU implementation via tensor broadcasting and dimension-we wise minimization, with infeasible transitions masked via +oc. For variable-sized state spaces, padding and masking ensure regular tensor shapes for parallel execution. A demonstrative DP example (A.2), together with the full formulations of the following two applications, CVRPSI) split (A.3) and DSIRP reinsertion (.4), is provided in the Appendix.
2.3
INSTANTIATION A: SPLIT DP ON A GIANT TOUR IN THE VEHICLE ROUTING PROBLEM WITH STOCHASTIC DEMAND.
Problem Motivation. In vehicle routing postprocessing, common task is to split "giant tour" [oi...., Gn.] into capacity-feasibl routes. Given demands Lh under scenario w and vehicle capacity Q, define state having served customers di d;. An action p i ends the previous route at p. starting a new one from p-1 to i. The departure depot is denoted by 0, and the destination depot by n+1. The DP explores all possible cut points p that define where start new route, and accumulates the minimal total travel cost for serving customers to i. (see the Figure in A.3 for better understanding).
Forward DP Recursion and Matrix Form.
The cost of serving subroute [opt1, 0,]
( = C0,0p11 +
Cor.aki1 + Co.ntll
Amp+1
which is feasible only E1
kkppt Cr I Q. We define masked transition entries as
(p,p,
if p and Ei Che Q,
(p, =
+X.
if < and Ek=p+1 Lo
Q (capacity violated),
t0c,
if 2 ( not applicable).

---

Published as a conference paper at ICLR 2026
Let (0) 0 and (i) be the optimal cost to reach state i. The forward-DPP update is then
(i)
min m (p) + (p, i)},
i= 1, 1=......
p<t
which is equivalently expressed as the masked min-plus reduction
(1) = (a)) @ (0:1-1,,,
(3)
where & denotes the (min, +) semiring product and (,, is the i-th column with all infeasible or undefined entries masked by +00 (see Appendix A.3 for a numerical toy example).
Two-dimensional DP parallelism on GPU (scenarios transitions)
Custemer
Customer 1-1
Sccnurtos
Castomeri
Scenarios
Custameri
Customee i+1
Somariss
Customeri-1
Figure 1: 2D DP parallelism on GPU (scenarios x predecessors). Each row corresponds to destination state i, each block within row represents scenario 4, and colored bars indicate feasible predecessors p For each (i,w) pair, all predecessors are expanded parallel to form (J*(P) + (p.i) : p i), followed by a column-wise min-reduction over that yields (i).
2D Parallelism on GPU. From equation 3 the computation at a fixed destination state factorizes over the Cartesian product Ω x {p p i). We therefore exploit 2D GPU parallelism across scenarios wa and predecessors < Each thread computes one pair (w.p) by loading (p, i) and the partial cost (p), forming (p) + (", i). A warp-/block-level min reduction across p then yields J (i) for that scenario. Launching such kernels for all scenarios in parallel computes the masked min-plus reduction of (- i) against J (0:i=1). 1). As illustrated in Figure] | this structure maps naturally to GPUs: scenarios wa are parallelized across columns, predecessors are reduced within blocks, and rows (states i) advance independently along the DP frontier.
2.4 INSTANTIATION B:
FORWARD INVENTORY REINSERTION DP IN DYNAMIC STOCHASTIC INVENTORY ROUTING PROBLEMS.
Problem Motivation. In the stochastic inventory routing problem, delivery schedules are often determined at an aggregate level and then refined through local reinsertion moves: customer that risks a stockout is reconsidered, and new visits are inserted into existing routes or additional trips are scheduled. The key challenge that reinsertion must balance two competing effects: (i) earlier deliveries create higher holding cost and may cause inefficiency in vehicle loading, (ii) postponing deliveries increases the probability of future stockouts under adverse demand realizations. The goal is to decide, for each customer, when to replenish under uncertain demand expected routing, holding, and stockout costs. DP provides a natural way to resolve this tradc-off, as it captures the temporal coupling of inventory states and demand uncertainty.
to minimize
The decision process follows a two-stage stochastic optimization framework. In the first stage, delivery and routing plan is established for Day 1. Specifically, the model determines which customers to replenish and how much to deliver, subject to vehicle capacity and routing constraints. These decisions made prior the realization of demand and are identical across all demand
4

---

d.w The recursion starts from Published as a conference paper at ICLR 2026
scenarios. Once the demand over the entire planning horizon (Days to H) is realized, second- stage decisions are made adaptively from Day onward. These include scenario -eeaaidddeeeneent routing and replenishment actions that respond the realized demand in each scenario (see the Figure ??
in A.3 for better understanding).
The objective is minimize the total expected cost, which consists of components: (1) First- stage costs including the Day routing cost and the supplier's inventory holding cost; and (2) Second-stage costs, which vary across scenarios and include: inventory holding costs and stock- out penalties at customers at the end of Day I; routing and delivery costs from Day 2 to H; and inventory holding and stock-out penalties from Day H.
2025 which The DP operator used to solve this problem follows the formulation in (Zhao et al.
focuses single customer over planning horizon I...., .,, H, starting with imitial inventory T and capacity U: The customer's 's demand is uncertain and modeled by finite set of scenarios Ω,
where d: denotes demand on day / under scenario For each customer, the decision at each day consists of:
whether customer is visited at time
the delivery quantity q which can only take two values: of {0, U. ((--!--1) that is, either
no delivery or replenishment up to full capacity.
and which vehicle route is chosen to accommodate this visit.
These decisions are scenario-independent, i.e., the same schedule applies across all w E Q while in- ventory evolution scenario-dependent The state variable the end-of-day inventory I.t updated as
1/ = max 1-11~- + -d Vt, w.
Here the DP systematically evaluates both replenishment options (no delivery vs full OU delivery), propagating inventory states forward in time and accumulating costs. See Appendix A.4 for nu- merical toy example.
Forward DP Recursion and Matrix Form. At each day t. customer either receives no delivery
(af = 0) or is replenished up to capacity (g) = U.- - The per-stage cost for scenario w consists of two components: (1) routing and detour costs associated with sending q/- denoted (g();; and (2) customer-side inventory holding and stock-out penalties h (I!""))) evaluated at the end-of- day inventory 1!
Let C7 (I) denote the minimum expected cumulative cost up to day for customer under sce- nario w, given that the day-t starting inventory T. By construction, (.) is a piecewise linear function of the inventory state. The forward recursion is
Of out
min of (0 (0.00:;:111))
(G14) + Fo(e)) + (5((,+11g)
where the inventory state evolves as 1.+1, max (0, 1- + of the initial inventory before day 1: CT = with 10 = IP.
To enable GPU-friendly computation, collapse the action space into state-to-state transition matrix:
(I,))
min OE (0.0, max -4 )=]
(F((g) +))) +oo no feasible leads from to J.
In practice, evaluating each entry A A may itself require enumerating finite set of candidate route options (e.g., alternative reinsertion choices), in which case
A (1, min A (1r). rek
Rows correspond to today's 's starting inventory I, columns to tomorrow's 's inventory J, and each entry stores the minimal cost of transitioning from to J.
5

---

Published as a conference paper at ICLR 2026
Define the state space
= 22 :0:: ;
and collect (J) over S as vector J E R!!!!, The forward recursion then becomes a masked min-plus matrix-vector product:
J (A 2 J = min ress
(1, +
(4)
JES
3D Parallelism GPU. From the matrix form, the computation stage factorizes the Cartesian product Ω x {(T + J)} x R. where R denotes the set of candidate route options for each transition. We therefore exploit 3D GPU parallelism across scenarios w, state transitions - J, and route options
R. Each thread computes tuple (w, I->J,r) by loading AI,
and the partial cost (1). forming (4) A, (rr). A warp-/block-level min reduction first performed across route options r, then across predecessor states I, yielding (J) for each scenario w. Launching such kernels for all w in parallel realizes the batched column-wiss min -plus updates the recursion, while also vectorizing over alternative delivery routes.
Three-dimensional DP parailelism (scemarios x transitions decisions trunsition) Stage t-1
Inventory Level
Delivery
Route
(M)
(0)
Scemaries
20
Transtioo (i(iv)
Delivery
Stage
Route
Inventory Level
Seemarãos
Figure 2: 3D parallelism on GPU (scenarios x transitions X route options). Each layer corresponds to stage f, with nodes representing end-of-day inventory levels T. Colored edges denote feasible transitions I + J under scenario-specific demands d= For each tuple (w, r)) threads evaluate the cost contribution (1) + A, rr), combining routing overhead with holding and stockout penalties. A two-level reduction (first across route options r. then across predecessor states I) yields (J) per scenario. The figure highlights how GPU parallelism spans scenarios, transitions, route options, turning the DP recursion into fully batched min plus update.
EXPERIMENTS
In this section, we present numerical experiments demonstrating the effectiveness and practicality of the proposed matrix-fomn DP framework. Section [3.] examines the role of large scenario sets in improving solution quality, and Section 3.2 shows the scalability of our approach. Sections 3.3
and verify that these findings also hold for our GPU-based implementation, while Section 3.5
evaluates its fcasibility in terms of computational resources. Additional comparisons with baseline methods are provided in Appendices/A.66 A.7 with Aatehdpixxx.... Appendix/A..Jnnalyinnn analyyiig the effect of the number of evaluated first-stage solutions.
3.1 SCALING THE SCENARIO SIZE IN STOCHASTIC PROGRAMMING.
The theoretical properties of empirical risk minimization (ERM) under mild regularity conditions establish that SAA solutions may suffer from bias with small sample sizes but converge consistently

---

Published as a conference paper at ICLR 2026
toward the true optimum the number of scenarios grows, with asymptotic (1///m) conver- gence rate (see Appendix |A.5 practice, we conduct experiments the DSIRP, representative setting where demand distributions are complex and cannot be adequately captured by simple parametric families. Figure reports the empirical behavior of SAA estimators the number of scenarios increases.
.a.endiiiA.Aeeee for a formal statement). To examine how these properties manifest in
Specifically, Figures Ba and 3b highlight the bias effect: with only a small number of scenarios, the estimated systematically underestimates the true expectation. As the sample size grows, the estimator increases and gradually stabilizes, consistent with the theoretical consistency guarantee. Figures 3c and 3d further illustrate the convergence behavior. As the number of scenarios increases from hundreds to tens of thousands, the SAA estimate approaches the true optimum, and its variance decreases at the predicted ((1//m) rate. The log-scaled plot shows that improvements continue to accrue even at large sample sizes, underscoring that "more scenarios" consistently yield better estimates rather than reaching premature plateau.
In general, these findings demonstrate that when the underlying distribution of uncertainty is com- plex or unknown, as common in real-world, data-driven settings, restricting the scenario set to only few hundred or a few thousand is insufficient. Substantially larger scenario sets are needed to reduce bias and improve solution quality. Our GPU-accelerated DP framework makes such large- scale, aaiaenriven stochastic evaluation computationally feasible in practice.
Bas saaben dww.adl Josibniee
Bia trde sxrred Aknzad Jauibaiee
Canagener wift marbe sceuaio Oloar (irCaneggaee mianeeeer of d'cceaiir Gng bakk
scdel
Figure 3: Empirical behavior SAA estimators in DSIRP. Top: hias under different demand distributions.
Bottom: convergence with increasing scenario size.
3.2 SCALABILITY WITH THE NUMBER SCENARIOS.
The above results confirm that larger scenario sets statistically necessary to reduce bias and achieve consistency in stochastic programming. We next examine whether such scaling is compu- tationally feasible. We evaluate the efficiency of GPU-accelerated DP operators against CPU baselines on two representative tasks: (i) the split operator CVRPSD and (ii) the reinsertion op- erator in DSIRP. All implementations were written in C+/CUDA and tested on a machine with an AMD Ryzen 7 9700X CPU (8 cores) and an NVIDIA RTX 2080Ti GPU with 11 GB memory. The CPU baselines include (i) single-threaded implementation and (ii) multi-threaded imple- mentation with 8 threads. The GPU implementations exploit the 2D/3D parallelism described in Section 2.
Left (CVRPSD split DP). As the number of scenarios increases from 10 to 105, the single-thread CPU runtime rises sharply and reaches minutes 10 scenarios, while the 8-thread baseline shows only moderate relief before saturating due to synchronization and memory- -annory- limits. The GPU curve grows nearly linearly and remains in the seconds range even 10 scenarios, yielding about 80x speedup over single-thread CPU and 20x over the 8-thread baseline at the largest setting.
Right (DSIRP reinsertion DP). The effect far more dramatic. At x 10 scenarios, the GPU implementation attains roughly 9.3 x 10t speedup versus the single-thread CPU and 2.26 10a versus the multi-thread CPU (see callouts in the figure). These gains stem from: (i) 3D parallelism (scenarios x inventory transitions route options), (ii) high arithmetic intensity with coalesced memory access, and (iii) warp-/block-leew reductions that keep the Bellman minima on-chip.
Across both problems, GPU parallelization shifts the computational frontier for scenario-basee based eval- uation: near-linear scaling the number scenarios with order-of- ofnnnnut to five-orders-off magnitude speedups (problem-dependdent) This throughput is what enables the very large, data-

---

Published as a conference paper at ICLR 2026
(a) CVRPSD split operator
(b) DSIRP reinsertion operator
Figure 4:
Runtime comparisons of CPU and GPU implementations across different scaling regimes. Left:
scaling behavior under 10 scenarios for CVRPSD. Right: large-scale evaluation for DSIRP.
driven scenario sets used in our SAA experiments, directly supporting the statistical benefits docu- mented in the previous subsection.
3.3 IMPACT OF SEARCHING SCENARIO SET SIZE ON DECISION QUALITY
We next examine how the number of evaluated scenarios influences the quality of the first-stage decision. Specifically, we solve the problem under different scenario counts, ranging from to 10* (i.e., 100, 1,000). For each scenario setting, the obtained first-stage solution is evaluated fixed large out-of-sample evaluation set 10 scenarios. Figure 5 reports the out-of-sample cost achieved by the best observed solution CVRPSD instances: x-n128 and x-n105.
92236
323044
327942
127400
327040
93238
133438
137850
227036
123634
07220
223750
72/0
45.00
10
Figure 5: Out-lf-sample performance of first- stage solutions obtained with varying observed scenario settings. Larger evaluation set yield more robust and lower-cost solutions.
Figure 6: Quality of the best solution obtained at each time under fixed time budget. GPU consistently achieves better decisions duc faster evaluation and thus larger effective search effort.
Results. When observing only few scenarios (e.g., a single scenario), the resulting first-stage so- lution is severely biased and performs poorly. Increasing the scale of available scenarios consis- tently improves robustness, with significant gains observed throughout Figure Such performance confirms the theoretical insight that larger sample sizes reduce estimation bias in sample-average approximation. These results demonstrate that our GPU-based framework, by enabling the evalua- tion of tens hundreds of thousands of scenarios within practical runtimes, leads to significantly more reliable first-stage solutions compared to CPU-based methods that are restricted to only a few thousand scenarios.
8

---

Published as a conference paper at ICLR 2026
3.4 DECISION QUALITY UNDER FIXED TIME BUDGETS
Finally, we compare the decision quality obtained under identical wall-clock time limits across the three implementations: CPU single-thread, CPU multi-thread, and GPU. Each method is given fixed runtime budget during which the split algorithm splits giant tours to obtain first-stage solu- tions. For fairness, all approaches are evaluated on the same problem instance with 10 available scenarios. Figure 6 reports the best penalized cost obtained within the allowed runtime.
Results. With small time budgets, methods return feasible but suboptimal solutions, yet GPU already provides noticeable advantage. As the time limit increases, the quality gap widens: GPU produces solutions that are consistently closer the true optimum, while CPU single-thread stag- nates and CPU multi-thread improves only modestly. This matches intuition: faster scenario eval- uation allows the GPU to explore many more candidate first-stage tours within the same runtime, thereby improving the probability of discovering high-quality solutions.
These results confirm that beyond scalability, GPU acceleration directly translates into superior de- cision quality under realistic time constraints, making it particularly valuable in operational settings where decisions must made quickly.
3.5 GPU-USAGE ANALYSIS
Monery thage Sale Spare Cand Moneyy
10 Naeee)
www. Neaerro 3,08 910611
Figure 8: Relationship between GPU utilization and the number of scenarios/customers in CVRPSD in- stances.
Figure 7: Relationship between GPU utilization and the number of scenarios/I,, DSIRP instances.
To assess the viability of our method with respect to GPU memory requirements, we conduct dedi- cated experiments for both the GPU-based CVRPSD and DSIRP solvers. Figures and summa- rize the resulting GPU-memory consumption. Results for both CVRPSD and DSIRP show that GPU memory usage grows strictly linearly with the number of scenarios and the state-space size. Even large configurations, such 10t scenarios for CVRPSD (3.98 GB) 100,000 scenarios for DSIRP (1.66 GB), can fit easily within common GPU limits, implying that an 11 GB GPU can support several million scenarios. This confirms that GPU memory is not practical bottleneck; scalability limited by computation time rather than memory capacity.
CONCLUSION
We presented GPU-based framework for executing second-stage dynamic programs in a scenario- batched and multidimensional manner, enabling full-fidelity stochastic combinatorial optimization scales that were previously impractical. By exploiting parallelism across scenarios, DP layers, and action or routing choices, our kernels can evaluate over million uncertainty realizations in a sin- gle pass. Experiments on a vectorized split operator for stochastic vehicle routing and an inventory reinsertion dynamic program show that the method scales nearly linearly with the number of scenar- ios and yields substantial speedups over multithreaded CPU baselines, leading directly to stronger first-stage decisions under the same computational budget. Although the framework inherits GPU memory limitations and currently focuses on forward dynamic programs with additive costs, it opens several promising directions. Future work includes handling constrained and risk-aware objectives,

---

Published as a conference paper at ICLR 2026
designing GPU-native abstractions for irregular dynamic programs, integrating learned components while preserving numerical reliability, and extending the primitives beyond the specific problems studied here. We view this work early step toward hardware-aware algorithms for large-scale stochastic discrete optimization.
ACKNOWLEDGMENTS
This work was supported by the Longgang District Special Funds for Science and Technology In- novation under Grant LGKCSDPT2023002.
REFERENCES
Cynthia Barnhart, Ellis L. Johnson, George I. Nemhauser, Martin WP Savelsbergh, and Pamela H Vance. Branch-and-price Column generation for solving huge integer programs. Operations research, 16-3);316-329, 1998.
Richard Bellman. On routing problem. Quarterly of applied mathematics, 16(1):87-90, 1958.
Richard Bellman. Dynamic programming treatment of the travelling salesman problem. Journal of the ACM (JACM), 91):611-3, 1962.
Dimitri Bertsekas. Dynamic programming and optimal control: Volume I, volume 4. Athena scien- tific, 2012.
Dimitri P Bertsekas and John N Tsitsiklis. Neuro-dynamic programming: overview. In Pro- ceedings of 1995 34th IEEE conference on decision and control, volume pp. 560-564. IEEE, 1995.
Dimitris Bertsimas and Ramazan Demir. An approximate dynamic programming approach to mul- tidimensional knapsack problems. Management Science, 48(4))550-565, 2002.
John Birge and Francois Louveaux. Introduction to stochastic programming. Springer, 1997.
Arun I. Bishop, John Z. Zhang, Swaminathan Gurumurthy, Kevin Tracy, and Zachary Manchester. Relu-qp: A gpu-accelerated quadratic programming solver for model-predictio control. In 2024 IEEE International Conference on Robotics and Automation (ICRA), pp. 13285-13292. IEEE, 2024.
Marco Antonio Boschetti, Vittorio Maniezzo, and Francesco Strappaveccia. Route relaxations gpu for vehicle routing problems. European Journal of Operational Research, 258(2):456-466, 2017.
Nadjib Brahimi and Tarik Aouam. Multi-item production routing problem with backordering: milp approach. International Journal of Production Research, 54(4):1076-109, 2016.
Robert L Carraway, Thomas L Morin, and Herbert Moskowitz. Generalized dynamic programming for stochastic combinatorial optimization. Operations Research, 37(5):819-88) 829, 1989.
Diego Cattaruzza, Nabil Absi, Dominique Feillet, and Thibaut Vidal. A memetic algorithm for the multi trip vehicle routing problem. European Journal of Operational Research, 236(3):833-848, 2014.
Geoffrey G Chu. Improving combinatorial optimization. University of Melbourne, Department of Computer Science and Software Engineering. 2011.
Leandro Callegari Coelho, Gilbert Laporte, and lean-François Cordeau. Dynamic and stochastic inventory-routing. CIRRELT Montreal, 2012.
Timothy Davis. Algorithm 1000: Suitesparse: Graphblas: Graph algorithms in the language of
sparse linear algebra. ACM Transactions on Mathematical Software (TOMS), 45(4):1-25, 2019.
Edsger W. Dijkstra. A note on two problems connexion with graphs. Numerische Mathematik, 1: 269-271, 1959. doi: 10. 1011/000011866600 386390.
10

---

Published as a conference paper at ICLR 2026
Joseph Farrington, Kezhi Li, Wai Keong Wong, and Martin Utley. Going faster to further: Gpu-accelerated value iteration and simulation for perishable inventory control using jax. arXiv preprint arXiv:2303.10672, 2023.
Dominique Feillet, Pierre Dejax, Michel Gendreau, and Cyrille Gueguen. An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems. Networks: An International Journal, 44(3):216-22) 2004.
Daniele Ferone, Paola Festa, Serena Fugaro, and Tommaso Pastore. A dynamic programming al- gorithm for solving the k-color shortest path problem. Optimization Letters, 1566)117731192, 2021.
Nicolai Fog Gade-Nielsen. Interior point methods on gpu with application to model predictive control. 2014.
Michael Held and Richard M Karp. The traveling-salesman problem and minimum spanning trees: Part ii Mathematical programming, 1(1):6-25, 1971.
Stefan Irnich and Guy Desaulniers. Shortest path problems with resource constraints. In Column generation, pp. 33-65. Springer, 2005.
Felix Liu, Albin Fredriksson, and Stefano Markidis. A gpu-accelerated interior point method for radiation therapy optimization. arXiv preprint arXiv:2405.03588, 2024.
Haihao Lu and Jinwen Yang. cupdlp. jl: A gpu implementation of restarted primal-dual hybrid gradient for linear programming in julia. arXiv preprint arXiv:2311.121800 2023.
Haihao Lu and Jinwen Yang. cupdlp- jl: A gpu implementation of restarted primal-dual hybrid
gradient for linear programming in julia. Operations Research, 2025.
Chryssi Malandraki and Robert B Dial. A restricted dynamic programming heuristic algorithm for the time dependent traveling salesman problem. European Journal of Operational Research, 90 (1):45-55: 1996.
Saeed Maleki, Madanlal Musuvathi, and Todd Mytkowicz. Parallelizing dynamic programming through rank convergence. ACM SIGPLAN Notices, 49(8):219-232, 2014.
Vangelis Th Paschos. Applications of combinatorial optimization, volume 3. John Wiley & Sons, 2014.
Warren B Powell. Approximate Dynamic Programming: Solving the curses of dimensionality, vol- ume 703. John Wiley & Sons, 2007.
Christian Prins. simple and effective evolutionary algorithm for the vehicle routing problem. Computers & operations research, 31(12):1985-2000 2004.
Giovanni Righini and Matteo Salani. New dynamic programming algorithms for the resource con- strained elementary shortest path problem. Networks; An International Journal, 51(3):155-170, 2008.
Michel Schubiger. Gpu acceleration of admm for large-scale convex optimization. Master's thesis, ETH Zurich, 2019.
Michel Schubiger, Goran Banjac. and John Lygeros. Gpu acceleration of admm for large-scale quadratic programming. Journal of Parallel and Distributed Computing. 144:55-67, October 2020. ISSN 0743-7315. doi: 10. 10.1016/j (0026/..pcc..202.055.022 URI. http:///dd wttp:/ooo...o/1.../ org/10. 1016// jpdc. 2020. 05.
Alexander Shapiro. Monte carlo sampling methods. Handbooks in operations research and man- agement science, 10:353-425, 2003.
Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczynski. Lectures on stochastic program- ming: modeling and theory. SIAM, 2021.

---

Published as a conference paper at ICLR 2026
Jeremy F Shapiro. Dynamic programming algorithms for the integer programming problem i: The integer programming problem viewed as a knapsack type problem. Operations Research, 16(1): 103-121, 1968.
Byung Duk Song, Kyungsu Park, and Jonghoe Kim. Persistent uav delivery logistics: Milp formu- lation and efficient heuristic. Computers & Industrial Engineering. 120:418 428, 2018.
JER Staddon. The dynamics of behavior: Review of sutton and barto: Reinforcement learning: An introduction Journal of the Experimental Analysis of Behavior, 113(2), 2020.
T. Vidal, T.G. Crainic, M. Gendreau, N. Lahrichi, and W. Rei. hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research, 60(3):61 1-64,, 2012.
Thibaut Vidal. Split algorithm in o (n) for the capacitated vehicle routing problem. Computers & Operations Research, 69:40-47, 2016.
Carl Yang, Aydın Buluç, and John D Owens. Graphblast: A high-performance linear algebra-based graph framework on the gpu. ACM Transactions on Mathematical Software (TOMS), 48(1):1-51, 2022.
J. Zhao, M. Poon, Z. Zhang, and R. Gu. Adaptive large neighborhood search for the time-dependent profitable dial-a-ride problem. Computers & Operations Research, 147:105938, 2022.
Jingyi Zhao, Claudia Archetti, Tuan Anh Pham, and Thibaut Vidal. Large neighborhood and hybrid genetic search for inventory routing problems. arXiv preprint arXiv:2506.03172) 2025.
A APPENDIX
A.1
RELATED WORKS
DP has long been a foundation of combinatorial optimization and stochastic control Carraway et al.
(1889] Chu 2011)
Paschos (2014), underpinning exact algorithms for the traveling salesman prob-
lmm Bellman 1962 Malandraki & Dial 1996), knapsack variants Bertsimas & Demir (2002). and
shottest-path computations Ferone et al. 2021 ) Righiii & Salani (2008). In vehicle and inventory
routing, DP- base split and labeling procedures are core subroutines in state-of-the-ae of-the-art heuristics and hybrid genetic/metabeurrsttess Vidal (2016); Cattaruzza et al.] (2014), enabling fast evaluation and repair of giant tours and complex neighborhoods. In this sense, DP forms the backbone of both exact and heuristic combinatorial solvers. MILP formulations offer another widely used modeling
approach for VRPs Brahimi & Aouam (2016); [Song et al. (2018), but solving such models using CPU-based solvers remains computationally intensive, especially for stochastic VRPs with large scenario sets.
Moreover, GPU acceleration has gained traction optimization and control, with notable devel- opments in GPU-ADMM SShuubieer et al. (220)), primal-dual methods for linear and quadratic programming Lu & Yangl (2025), and GaU-based interior-pone point algorithms Gade-Nielsen (2014). These works demonstrate that careful kernel design and memory organization can deliver substan- tial speedups relative to multithreaded CPU implementations. Researchers have also explored GPU acceleration for routing problems, either through algorithms specifically tailored GPU architec- tures or by leveraging high-level GPU libraries. For instance, BBsshhtti et al. (2017) accelerates q-route-type DP relaxations for VRPs GPUs, while MMalkii al. (2014) introduces a parallel linear linear-tropical DP algorithm that circumvents cross- stage dependencies via rank convergence. On
the library side, Davis (2019); Yang al. (2022) implement GraphBLAS for sparse linear algebra on GPUs, and Farrington et al. (2023) uses the GPU-enabled JAX framework to accelerate DP-based stochastic control and routing computations.
Despite these advances, current approaches lack structural generality. GPU-oriented libraries such as JAX require reformulating problems into linear-algebra primitives to match the library's 's abstraction, while algorithm-specifi GPU designs rely on bespoke kernels and custom data layouts, limiting portability across problem classes. To the best of our knowledge, no unified framework currently exists for accelerating based stochastic combinatorial optimization GPUs. A general- -prrpose framework would make GPU acceleration reusable across a broad range of Da-based problems, eliminating the need to redesign complex kernels for each new application.
12

---

Published as a conference paper at ICLR 2026
A.2 ILLUSTRATIVE EXAMPLE FOR THE MIN-PLUS LATRIX-VECTOR BELLMAN UPDATE.
To illustrate the min-plus matrix-vector Bellman update, consider a dynamic programming recursion over three states in stage and two states in stage (+
Let the cost-to-go vector at stage /. under scenario w be:
J =
3
and the transition cost matrix from stage t t + be:
22
5
AV =
+00
3
0
where (i denotes the cost of transitioning from state at stage i. to state j at stage + 1. The value +o0 denotes infeasible transition.
Then, the Bellman update in the (min, +) semiring becomes:
J = (At) @ J
which is computed entry-wise as:
J (1) min (17 J (1)) min 1 1 1,1,3+++33 m,,22, =2
J4 (2) min (4,
min + t00 +1 0 3] min +00,
Thus, the updated cost-to-go vector at stage t is:
2 3
J
This computation reflects forward shortest-path propagation layered graph, where the cost of reaching each state in the next stage determined by minimizing over all incoming transitions from the previous stage.
Algorithm High-level Scenario- Batched GPU Forward DP
Input: Scenarios 12, horizon T, state sets {St), transition-cost generator Ay, initial cost-to-go J
Output: Final cost-to-go vectors {J Jwen, t=1..7
1: Preprocess all scenario tensors (pad, mask infeasible transitions, move to device)
2: for t = 1 to T 1 do
3: Launch GPU kernel over the chosen parallel axes:
scenarios
state transitions (s ++8) 8') E St x S
optional actions Ac(s) (if present)
4:
Each thread:
load (s)
load transition cost s';
compute partial cost e= J A (s,s(a))
5: Within warp/block:
reduce min over actions a (if applicable) reduce min predecessor states write (s')
6: end for
7: return All Jr
A.3 INSTANTIATION A: SPLIT DP ON GIANT TOUR IN THE VEHICLE ROUTING PROBLEM WITH STOCHASTIC DEMAND.
Problem Statement. The stochastic programming community has extensively studied scenario- based formulations, where uncertainty is modeled by finite set of realizations. However, scenario- based evaluation quickly becomes computationally prohibitive on CPUs, where even tens of thou- sands of scenarios overwhelm multi-threaded implementations. In work, adopt this
13

---

Published as a conference paper at ICLR 2026
scenario-based modeling framework, in which customer demands are represented by sampled alizations. This approach naturally accommodates correlated demand structures and supports data- driven modeling when historical records are available. This rrenario-bsse approach falls natu- rally into the two-stage stochastic programming paradigm, whose general form is: minzex fi (x) + Egg (fa(a, where denotes the first-stage routing decisions and the corresponding cost fi (x)) is a random vector representing a realization of customer demands, and f2 {) denotes the second- stage recourse cost under each scenario Ç
In our two-stage stochastic optimization problem, the first stage determines the visiting sequence of customers, commonly referred as a giant tour in the context of genetic algorithms (Vidal et al. 2012). This representation encodes solution as a permutation of customers, from which feasible ve- hicle routes can recovered via a split operator under fixed vehicle capacity and scenario-depeenent customer demands. We assume full demand revelation prior to the second stage, enabling the plan to adapt to the realized scenario. Once demands are known, the giant tour is split into feasible routes such that the demand on each route does not exceed vehicle capacity. Thus, given a fixed giant tour (i.e., the visiting sequence), the second-stage evaluation computationally simple: it only requires splitting the tour according realized demands. The objective is to determine first-stage giant tour that minimizes the expected total travel cost across all scenarios.
First stagy
Chromosome
3
Second
Secnario
Scenasio
Scenario
119
4
2
2
4
3
Figure 9: Example of splitting giant tour into feasible routes under a given demand scenario.
The extensive model is:
min 2,
(8ig
- Sy) 2ijA + Σ 1
(5)
Ω
WER
REA
(KJ)
AEK
s.t.
ΣΣ
Titk
Vj E W (0}, Vo
(6)
KENSEN
Σ
19k =1
V W (0), Ae
(7)
JEN
Σ
Σ
E (0), V E K. Vo
(8)
jikk
iEN
VEN
Ak K, AM
Oojk
1#0
170
Σ
S
K
(10)
1
14

---

Published as a conference paper at ICLR 2026
ΣΣ JEN
a Di - C S 1
Ak K, PM
(11)
4/ Ujk +15 n(1 xz(zk)
# J. 1,j 2 1, Ak, Ae
(12)
Σ aj + Yıj= =1
(13)
Σ Fij + Yij =1
W22
(14)
26 Σ
Fijk
vi> 0, m
(15)
2ij 2 Σ
10k + Σ
Oikk - 1- - Sij
0, 0, m4
(16)
k
b
4j U1 U, -1
vi,J
(17)
2;2 U, - -U, +1
vi.j
(18)
U- - U, +1 n(1 - zzz)
(19)
ΣΣ Yg =1
(00)
ziju E (0.1) Zij E {0.1)
1i サ J. Ak, Ye サi Vi, k,
(21)
(22)
1 1 u Sn
(23)
12 120
Ak, w
(24)
Suppose the first-stage giant tour is Figure
with vehicle capacity Q = 17, as illustrated in
(i) In Scenario 1, the realized demands are [14, 15. 8, 1. The second-stage routing may then split the tour into three feasible routes: (1), (2), and (4.3, 5).
(ii) In Scenario m, the realized demands are [8, 1, 5, 7, 2] with three vehicles available, The second-stage routing may then split the tour into two feasible routes: (1) and (2, 4, 5,5.
is well understood that searching robust first-stage policies benefits from incorporating a large number of demand scenarios to capture the full breadth of uncertainty (Shapiro (2003), see the details in next section. Theoretical results show that a sufficiently large scenario set is critical for achieving accurate and stable estimations. However, in practice, evaluating even tens of thou- sands of scenarios can be prohibitively expensive, especially for combinatorial problems such as the CVRPSD, since each scenario requires solving non-trivial routing evaluation (typically through dynamic programming algorithm with time complexity O(nP), where n is the number of customers and P is the number of possible transitions) rather than a simple function call.
Pseudo-code for 2D Kernel. In this section, we describe the GPU kernel designed to generate data splits for each scenario in Algorithm 3 In Algorithm 2 we launch this kernel with a grid of thread blocks to materialize the output matrix Splits for all scenarios are produced concurrently, which enables simultaneous split generation across all scenarios.
Illustrative Matrix Form of DP. To clarify the structure of the transition cost matrix A in the split problem, consider a simple instance with a giant tour [01, 022 0z), and realized demands under scenario w given by: ( Lo to (2, 4), with vehicle capacity Q 5.
Let (1) denote the minimum cost to serve customers dn through or, We construct transition cost matrix A e 3333 where the (p. i) entry represents the cost of serving customers 0p+1 to o; in one route, if the cumulative demand is within capacity.
15

---

Published as a conference paper at ICLR 2026
Algorithm 2 GPU Splitting Algorithm
Input: Scenarios vi E !
Input: Global settings Q, C. for vehicle capacity and travel cost matrix.
Output: Running costs V
1: Initialize V E R QI to all zero matrix
2: Instantiate multiple 2D kernels parallel for VW.
3: V = [Vwo, Vwi [wa]]]
4: return
Algorithm 2D Kernel for Splitting
Input: Scenario Index w
Input: Scenario-specific demand q"
Input: Global settings Q, C' for vehicle capacity and travel cost matrix.
Output: Running costs V
I: Initialize V- to all zero vectors
2: Initialize Ç to an empty queue.
3: for each customer do
4: Update potential V/
5: < n then
6: if i dominates Shack then
7:
while 5"d *0 AND Sback dominates from right do
8:
Pop q" from back
9:
end while
Push i X from back
11:
end if
while | > AND Airon better than (w
fnext_front do
13: Pop Ç" from front
end while
15: end if
end for
17:
return V'
Assume the travel cost matrix is:
0
1
37
1
[Co,p]
2
3
2
0
where node ( is the depot and node 4 is the return depot. Then we compute:
A 1): route [o1]. demand = < 5 feasible. Cost = Co,1 + c1,4
1
(0, 2): route [01, 02), demand 5 feasible. Cost =1+1+2 4.
(0,33)) route [3)) demand = 9 infeasible. Cost +00.
A"(1,2 route [02]. demand = feasible. Cost= c,,2 + 02,4 = 2 + 4.
(1,33) route [02, 03), demand infeasible, Cost +00.
(2,33) route [og]. demand 4 feasible. Cost CD,3 + C3,4 3 4.
Thus, the masked transition matrix becomes:
x
+oo
to0
A
X
x
the forward DP proceeds column by column from the initialization J²(() = 0:
J = (0)
16

---

Published as a conference paper at ICLR 2026
(since only the predecessor p = 0 is admissible when 1). Then
min P00
=
+ =
and
(3) min
+ (2 (F0 min +4+ 4 too, 0 + +00
Thus the running costs after this pass are
w
= (0, 4, 4,
This is exactly the left-to-right scheme you described: each new (J) is the minimum over the previous prefixes, i.e., J [ = min (r) + (r:jj)) so the recursion naturally continues as
j j+j++ on larger instances.
A.4 EVALUATION ON INSTANTIATION B: FORWARD INVENTORY REINSERTION DP IN DYNAMIC STOCHASTIC INVENTORY ROUTING PROBLEMS.
Problem Statement. We consider a two-stage stochastic version of the Inventory Routing Problem (IRP) under an Order-Up-To (OU) inventory policy. We define the model on complete directed graph G (N, where N {0, n + 1} A includes the supplier node 0, the destination depot + 1, and the set of customers .. The are set A represents all possible directed connections between nodes. set of vehicles IC is available for deliveries, with KC] = K. The planning horizon spans finite set of days T {1,.... ......,,) and T' = {2,... ..... H/} for the second stage. To model uncertainty in demand, we consider finite set of scenarios St where each scenario w Ω occurs with probability pd. Each vehicle has capacity limit Q, and each customer capacity U.. Holding costs h2 are incurred per unit of inventory stored at node i, and stock- out penalty incurred at a rate of phş per unit of unmet demand at customer i, where p >> 1. The customer demand at customer i day under scenario is denoted by d.w. The cost of traveling from node i to node j is denoted Gjj for each arc (i,j) E A. Finally, 7 represents the initial inventory level at node at the beginning of the planning horizon. To simplify replenishment decisions and reflect common logistics practices, we adopt an OU policy. Under this policy, each customer is either replenished up its full capacity U,, or replenished at all on given day. This f.w in the second stage), which equal if customer is is modeled using binary variables zद (or 2; replenished on day f, and 0 otherwise.
N' has an inventory
The first-stage of the 2SIRP involves determining the vehicle routing and delivery quantities on Day 1, before the actual demand realizations are observed. Given the first-stage decision and realized demand d= = {dd hEN.ETT under scenario w, the second-stage cost function Q (x)) can be calculated ho +kEE E(sjee GyyVej k.1 E Q²)). The second-stage constraints under the OU policy govern the evolution of inventory, replenishment quantities, and routing decisions from Day 2 to Day H under each demand scenario w E S2. Constraints enforce the OU policy: a customer either receives a shipment that raises its inventory up to the full capacity U; (i.e., q k.t.w
U, - U ---,",,, receives nothing. Constraints also need to ensure cother classical IRP constraints, such as the total quantity delivered by any vehicle i on day t does not exceed its capacity Q.
Pseudo-code for 3D Kernel. To efficiently adapt the CPU-based DP algorithm for GPUs. we design two-part framework. Specifically, Algorithm [ selects a GPU-feasible batch size, solves scenarios in batches via DP kernel, and aggregates results, reducing the batch size and retrying upon out-of-memory Algorithm executes a batched backward dynamic program over the hori- zon-parallel across scenarios and states-comparty recording the minimizer and transition, and then backtracking from each scenario's 's initial inventory to recover daily decisions, quantities, and total cost.
"no delivery" versus "deliver to capacity,"
Illustrative Matrix Form of DP. Consider single customer with U; 2, horizon 1, 2, and
single scenario w. The initial inventory is 7 daily demand is d:."" 1, the transportation cost
17

---

Published as a conference paper at ICLR 2026
Algorithm 4 GPU-Based DP with Adaptive Batching
Input: params, client_id, wal.sseeaarion iasararoos_dattt data. initial batch size Bu, device
Output: List scenario results (cost, decision flags, quantities)
1: B < ADIUSTBBTTCHSIIEE( Bn. device, params)
based on free GPU memory
2: results + []: start - 0; N + #scenarios
3: while start < NV do
4: batch < slice(all www.allsseaarrrreeeeeee is data, start:start+
5: (costs, flags, qt.) SOLVEBATCHDD((atch) params, device)
6:
per-seenrro tuples to results
7: start start +B
8: CLEARCACHEO
9: end while
10: return results
Algorithm 5 SOLVEBATCHDP
1: function BOLLEEBATCHDPPatt
params, device)
2: Preprocess tensors (time-major, contiguous, to device); set horizon T, state grid S Initialize DP arrays C, D, P: set C[T] for t = T 1........0 do
4:
(parallel over scenarios/states on GPU)
5:
Compute no-delivery cost and next state
6:
Compute deliver-to-me to- cost (fixed + capacity + holding) and next state
7: Cit] min(-); D argmin decision; P(t] next-state index
8: end for
9: return BACKTRACK(C, D, start inventories)
costs, flags, quantities
10:
end function
is F((g) q. and the end-of-day inventory cost is
2,
M(I
I /= 1,
5, ==0. O.
The state space each day is S {0, 1,)
Day 1. The forward transition matrix A! A!, J) (rows I, columns J) is
5
3
+00
1
2
to0
(g U. - 17, J max(0, 9 -1}
+oo
1
to0
The initial cost vector encodes TQ 1:
+
J
The column-wise min-plus update gives
min 0+5, 00) min 0+2, oc+1} min(co, 0o, 20) +o0
JR(J) min (1((1) re((.1.2)
A! J
Day 2. Parameters the same, A? A Updating once more:
min(5+5, 2+5, oo
7
(I, (IJ),
J = m5+3 2+2, c=+1)
=
4 +oo
min 0o, oo)
With two-day horizon, the terminal cost is min, (²(J) 4.
18

---

Published as a conference paper at ICLR 2026
(replenish full) is optimal. leading with cost 1+ 1 = Day 2: again from I = 1. delivering = is optimal, adding another 2. Total cost 2 2 4, which matches Jp. If Day 1 skips delivery, I drops 0 (cost 5), and even if Day delivers 2, the total becomes 8, which is suboptimal.
Policy insight. Day 1: from I 1, delivering
Notes. (i) The third column of A is always +oC because with = the end-of-day inventory cannot exceed 1. (ii) In general, if evaluating A(I,J) requires enumerating multiple route options E K. then A(I,)) min, A(1, J;r), and GPU parallelism naturally extends three dimensions (w, I-vJ, r) with reductions first over and then over I.
A.5 MONTE CARLO METHOD PROPOSITION.
Empirical Risk Minimization (ERM) method analogous the Monte Carlo method for estimating a population mean via sample averages and fits naturally within the ERM framework for stochastic programming. Formally, we distinguish between:
True problem:
(P)
=" [If O)
ZEX
Sample-average problem with m scenarios:
1 xE
(Pm) =m = min f(x
As the sample size m increases, the optimal solution v. of the sample problem converges to the true optimal solution z", and the optimal value 2,, approaches 2². The following result summarizes the fundamental properties of the ERM method under mild regularity conditions (cf. (Shapiro) 2003). ).
Proposition 1. Let (E ..'. be ii.d. samples of Denote by
m E
(P) min (I)(,r,)) () zE
(Pm. =
min xEX m.
(x,
Y=1
the true and sample average problems, respectively, with optimal solutions x and x. Then:
(Bias)
Bi 5) S =²
(25)
(Consistency)
BI E) ah x as m + 20,
(26)
(Probabilistic Convergence)
lim Pr B[f 5) = + ēmm} 2 a, ém 40, (27)
m=0
(Rate of Convergence)
Vm (=m
2²) Sn o0,000 oz())))
(28)
A.6 GPU BASED SPLIT VS. BASELINE MILP ON CVRPSD
We the vate-of-the-ar MILP solver Gurobi solve the CVRPSD extensive form under varying numbers of scenarios and customers, and compare its performance with our GPU-hased splitting algorithm. The results presented in Table show clear scalability gap: while GPUSPLIT'' time grows roughly linearly with the number of scenarios, the CPU-based MILP quickly becomes impractical, hitting time limits and eventually exhausting memory. In addition, for larger scenario sets, GPUSPLIT remains tractable and often achieves better objective values. For example, on the 10-customer dataset with over 1,000 scenarios, neither method proves optimality, but GPUSPLIT consistently finds higher-quality solutions than the CPU-based solver.
A.7 GPU REINSERTION VS. BASELINE ALGORITHM ON DSIRP
Regarding DSIRP, we provide comparison with the state-of the-aar method in Coelho et al. (0012), which solves DSIRP under single scenario. As shown in Table I across instances of various sizes, our approach achieves 22%%-6377 lower objective values (5.09% on average), while keeping run- times in the same order of magnitude and even running faster the largest class. This demonstrates that explicitly handling multiple demand scenarios can yield consistently higher-quality DSIRP so- lutions without incurring prohibitive additional computation compared to Coelho al. (0012).
19

---

Published as a conference paper at ICLR 2026
Table 1:
Comparison of GPUSPLIT and the MILP solver on CVRPSD instances with different
numbers of scenarios, reporting runtimes, objective values, and failure cases (time limits and 0OM).
#Scenario
100
500
1000
2000
10000
Time (s.)
GPUSPLIT
13.76
28.92
52.64
95,42
420.11
MILP
26.7
587.33
1708.51
3600
OOM
6 Customer
Obj. Value
GPUSPLIT
3399.59
3392.72
3391.86
3391.43
3391.17
MILP
3399.59
3392.72
3391.86
10541.3
Time (s.)
GPUSPLIT
485.03
3215.46
3600
3600
3600
MILP
3600
3600
3600
3600
OOM
10 Customer
Obj. Value
GPUSPLIT
3797.37
3804.89
3800.62
3802.95
3955.16
MILP 3797.37
3804.89
3913.01
5336.41
Table 2:
Performance comparison between the baseline DSIRP algorithm and multi-scenario method in terms objective value, runtime, and percentage improvement.
Instance
Baseline Obj.
Baseline Time (s)
Ours Obj.
Ours Time (s)
Gain (%)
small (n < 50)
9131.50
67.10
8703.20
102.30
4,68
medium (50 < n S 100)
30137.80
888.30
28219.00
1029.80
6.37
large (n 100)
60051.40 1.00
9248.90
57523.30
7200.00
4.21
Average
33106.87
3401.40
31481.83
2777.40
5.09
A.8 RESULTS ON SOLUTION QUALITY NUMBER EVALUATED CANDIDATE FIRST-STAGE SOLUTION
Regarding the number of candidates evaluated, we have conducted experiments to address this ques- tion, and present the numerical results in Table The quality of first-stage decisions number of candidates is an important consideration. We have presented partial results in Figure 6: the heuristic runs longer, discovers better first-stage solutions, which be interpreted as evaluating a larger and more diverse set of candidates. To make this point more explicit, we present an additional of experiments on the DSIRP instance, which directly reports the evaluation of different first-stage solutions and their corresponding performance under large scenario sets.
varied
Table Relationship between explored candidates and objective improvement in the DSIRP dataset.
Time (s.)
GPU Calls
# Candidates
Obj-
Gain
120
1,380
92
77.898,383
240
2,715
181
77,840,004
0.07% J
360
4,065
271
77.804,713
0.12% |
480
5,415
361
77,741,966
0.20% _
600
6,750
450
77,741,966
0.20% 1
The candidate first-stage decisions are generated by a genetic algorithm and evaluated via GPU, where each GPU call solves 500-scenario DP subproblem in parallel. Based on the result in T on the 15-customer, 500-scenario instance, this results in 92-450 evaluated candidates de- pending on the time budget, during which the framework consistently finds equal better solutions as more candidates are considered. The reported results clearly show how solution quality evolves with the number of evaluated candidates, thereby directly illustrating the behavior and viability of our approach under varying candidate counts.
We conducted the same set of experiments for CVRPSD (Table 4) and observe the same behavior as in DSIRP: as the time budget increases and more first-stage candidates are explored, the framework consistently returns equal better solutions, with no degradation in solution quality. The exper-
20

---

Published as a conference paper at ICLR 2026
Table 4: dataset.
Relationship between explored candidates and objective improvement in the CVRPSD
Time (s.)
GPU Calls
# Candidates
Obj.
Gain
120
304,740
304,740
3,893.32
240
611,841
611,841
3,811
2.11%
360
914,889
914,889
3,800.62
2
480
1,216,803
1,216,803
3,800.62
2.38%
600
1,524,407
1,524,407
3,800.62
2.38%
imental results clearly track how the objective improves with the number evaluated candidates, further confirming the robustness and viability of our candidate- baaddddte-baeed GPU evaluation scheme.
21


Reviews:
[{"review_id": "oahnQZk0IE", "reviewer": "Reviewer_7hTf", "summary": "This paper proposes a novel approach to stochastic programming problems (with DP representable recourse) that is highly parallelizable on GPUs.   Computationally, the authors demonstrate the efficacy of their approach on large-scale inventory and routing problems.", "questions": "- Is there any way to extend this to problems that are not representable as DPs?\n- How many candidate first-stage solutions are used?\n- How does the performance vary with the number of first-stage decisions?\n- The word 'training' is somewhat misleading, as it implies a notion of training that is not present.  In stochastic programming, it is essentially solving over different-sized scenario sets.  \n- I believe the authors may be using the ICLR 2025 template.\n- The code is not anonymized.", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- **Motivation**: Overall, the evaluation of candidate solutions over many scenarios in stochastic programming is computationally challenging, and the authors propose a unique new perspective that may have relatively strong applicability within the framework proposed, but also in many other settings.  \n- **Efficiency**: The method is very efficient, especially for problems with large-scale scenario sets and complex second-stage problems.", "weaknesses": "- **Scope & Generality**: One major drawback of this approach is the requirement of being able to represent the second-stage problem on a GPU through a DP.  While this does capture some stochastic programming problems, it may not be as applicable as other classical/ML-based methods.  Another weakness of the process is the requirement to implement a DP formulation for each new problem.  Compared to existing stochastic programming methods that formulate second-stage optimization problems, the barrier to implementation would be significantly higher. That said, there are certainly use cases where investment in implementation would be worthwhile.   \n- **Baselines**: This paper explores relatively large-scale benchmarks and demonstrates the efficiency of their approach.  However, their approach essentially acts as an evaluation on a set of candidate first-stage decisions, rather than truly optimizing over the entire decision space (what most methods do in stochastic programming). For that reason, it isn't easy to assess the quality of the solution in this case.  Comparing against established exact methods (integer L-shaped, Benders' dual decomposition) or heuristics (progressive hedging via MPI-SPPY) would help assess decision quality.  It may only be possible to get these methods working on smaller-scale instances, but it would still be highly informative.  \n- **Clarity**: Some details in this work are a bit unclear.  For example, based on my understanding of the paper, they obtain a set of candidate first-stage decisions and then use GPU parallelization to evaluate them; however, it is unclear how these decisions are received and how many are evaluated.  Moreover, there is no exploration into the decision quality over a varied number of first-stage candidates, which is necessary to assess the viability of this approach.  \n- **ICLR Fit**: While this paper explores GPUs, which are widely used in ML, it doesn't actually do any ML, only exact evaluations through parallel DP operations.  For that reason, there may be more suitable venues, such as CPAIOR and INFORMS Journals, etc.  That said, there is a growing interest at the intersection of stochastic programming and ML, so there is some relevance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "WVAJsdzqnO", "reviewer": "Reviewer_df1c", "summary": "This paper addresses the scalability bottleneck of Dynamic Programming (DP) in large-scale Stochastic Programming (SP) problems. A general-purpose GPU framework is proposed, which reformulates forward DP recursions as \"batched min-plus matrix-vector products\" over layered Directed Acyclic Graphs (DAGs). This approach allows the computation to be mapped directly to GPU kernels, enabling massively parallel computation. Validated on two classic SP applications, the framework demonstrates the capability to process large-scale scenarios, achieving a speedup compared to traditional multi-threaded CPU baselines.", "questions": "The acceleration at the operator level is significant. To evaluate the practical impact of this work, could the authors integrate this GPU operator into a complete solving framework and compare its end-to-end performance (in terms of total runtime vs. solution quality) against a SOTA solver on the same CVRPSD/DSIRP instances?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. This paper addresses a critical gap by extending GPU acceleration from its mature application in continuous optimization solvers to the computationally intensive domain of SP.\n2. It introduces a general-purpose GPU framework that reformulates broad classes of forward DP recursions into batched min-plus matrix-vector products. \n3.  The framework achieves a speedup compared to traditional multi-threaded CPU baselines.", "weaknesses": "1. The experimental comparisons in this paper are conducted at the operator level rather than the problem level. The core performance graph compares the authors' GPU DP kernel against their own CPU DP implementation. This constitutes a relatively weak baseline. While it demonstrates the operator's throughput, it fails to prove that the end-to-end performance is superior to a SOTA solver when this operator is integrated into a complete solution framework.\n2. The paper lacks external benchmarking for solution quality. Although Figures 5 and 6 effectively illustrate internal trends, the study never compares the final solution quality found by its method against publicly available \"best-known solutions\" in the literature. We learn that the operator is fast, but it remains unclear whether the solutions it finds are competitive.\n3. The paper lacks a dedicated \"Related Work\" section, making it difficult for reviewers to clearly position and contrast this work's contributions against other SOTA methods in the field.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "if9Q0tnJDn", "reviewer": "Reviewer_8Uwb", "summary": "This paper accelerates finite-horizon forward dynamic programs by formulating each Bellman update as a batched min-plus matrix–vector operation and is engineered to run in parallel on GPUs. Two applications are demonstrated: (1) a VRP split problem with stochastic demands (choosing route cuts along a fixed giant tour), exploiting 2D parallelism across scenarios and transitions, and (2) a dynamic stochastic inventory routing reinsertion operator, exploiting 3D parallelism across scenarios, inventory transitions, and route options. The results show higher scenario throughput and improved expected-cost decisions under the same wall-clock budget.", "questions": "1. You present the min-plus matrix–vector product as a core contribution. Please clarify how this differs from prior work such as [4] (DP as a tropical matvec) and existing semiring frameworks like [5] (references refer to the ones mentioned in the \"Weaknesses\" box\n\n2. How does this work differ from earlier GPU-based dynamic programming methods?\n(a) VRP: comparison to [1], who used GPUs for DP-based q-route and ng-route relaxations.\n(b) Inventory: comparison to [3], who used GPU-based value iteration (parallel Bellman updates) for large state spaces.\n\n3. The DSIRP speedup (9.3 × 104) looks extremely high. please clarify how the CPU baseline was implemented, in particular:\n(a) What algorithm was it implementing?\n(b) Why didn’t you include a CPU implementation of your own min–plus algorithm for a fair comparison?\n(c) Could you also compare with a semiring library baseline (e.g., [5]) to separate hardware from algorithmic effects?\n\n4. The framework seems to require, in the dense worst case, batched transitions of size O(|Omega|m2) with |Ω| = 106.\n(a) How was this handled on an 11GB GPU? Please report the per-stage state sizes mt used in CVRPSD and DSIRP.\n(b) Please include a memory-complexity analysis showing the trade-off between |Omega| and m. If large scenario counts are only possible for very small m (e.g., m < 50–60), how does that affect the claimed generality of the framework?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "1. The paper neatly reformulates forward DP over layered DAGs as masked min–plus matrix–vector products, making the recursion uniform, parallel, and easy to adapt.\n\n2. The use of masking with warp- and block-level reductions, and structured 2D/3D parallelism is technically solid and makes effective use of GPU hardware.\n\n3. Experiments demonstrates how large-scale scenario batching achieves high GPU throughput, clearly linking faster computation to more stable estimates and better decisions under fixed time budgets.\n\n4. The VRPSD and DSIRP examples show that the approach applies to different DP structures, and the toy examples and algorithmic sketches improve clarity.", "weaknesses": "1. Overstated novelty and missing literature. The paper treats the “batched min–plus matrix–vector” view as a key contribution, but this idea is well known from tropical/semiring DP and existing GPU semiring frameworks [4, 2, 5]. It should also acknowledge prior GPU-DP work in [1] and [3]. In my mind, the true novelty lies in the scenario-batched kernels and problem-specific engineering, not in the algebraic formulation itself.\n\n2. The GPU is compared to a different CPU algorithm, mixing hardware and algorithmic effects. A multi-threaded CPU version of the same method would make the comparison more fair.\n\n3. The paper notes GPU memory limits but gives no quantitative analysis. With 106 scenarios on an 11 GB GPU, state grids per stage must be very small if transitions are stored densely. This suggests that large scenario batching trades off against state-space size, a limitation that should be discussed.\n\n4. The VRPSD and DSIRP examples are solid but narrow; extending to more complex DP forms (non-layered, constrained, or\nrisk-sensitive) remains unclear. I feel that the work is better seen as a systems contribution for forward DP rather than a general GPU-DP framework.\n\nReferences\n\n[1] M. A. Boschetti, V. Maniezzo, and F. Strappaveccia. Route relaxations on GPU for vehicle routing problems. European Journal of Operational Research, 258(2):456–466, 2017.\n\n[2] T. A. Davis. Algorithm 1000: Suitesparse: Graphblas: Graph algorithms in the language of sparse linear algebra. ACM Transactions on Mathematical Software, 45(4):44:1–44:25, 2019.\n\n[3] J. Farrington, K. Li, W. K. Wong, and M. Utley. Going faster to see further: GPU-accelerated value iteration and simulation for perishable inventory control using JAX. arXiv preprint arXiv:2303.10672, 2023.\n\n[4] S. Maleki, M. Musuvathi, and T. Mytkowicz. Parallelizing dynamic programming through rank convergence. In Proceedings of the 19th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (PPoPP’14), pages 219–232. ACM, 2014.\n\n[5] C. Yang, A. Buluc, and J. D. Owens. Graphblast: A high-performance linear algebra-based graph framework on the GPU. ACM Transactions on Mathematical Software, 2022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uVUcgxclwp", "reviewer": "Reviewer_fYpR", "summary": "The paper presents a novel and impactful GPU-based reformulation of dynamic programming (DP) for large-scale stochastic combinatorial optimization. By expressing DP recursions as batched min–plus matrix–vector products over layered DAGs, the authors unlock massive scenario-level and transition-level parallelism. The approach demonstrates impressive speedups and enables practical evaluation of millions of uncertainty realizations, leading to improved SAA solution quality under realistic time budgets. However, there are several minor issues and open questions that need clarification or further validation. If the authors can address these points during the rebuttal, I would consider raising my score. The contribution is solid, and the direction is promising.", "questions": "1. The authors talk about the limitation of the approach. A related question: will the approach still be applicable if the action or the state is continuous? \n2. Could you explain what will happen in the edge case, like if there’s some irregular state spaces, how does the algorithm perform.\n3. From the 2D and 3D example, it seems like the algorithm design is a bit different in the appendix. Is it possible to write a unified algorithm and somehow put it in the main body instead of the appendix? Since the specific algorithm steps are very important for the paper. \n4. The paper uses DAG several times, and I suppose it refer to Directed Acyclic Graph, but I don’t think it is formally defined (Please ignore this if I am wrong). Also, some intuitions about why DAG matters here can make the paper easy to follow.\n5. Some terminology inconsistency. E.g., two-dimensional vs 2D.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "1. The idea is novel: By formulating stochastic programming (SP) as a specially DP, the paper proposes a novel approach which leverages GPU to solve SP. While the traditional sampling-based SP only includes thousands of scenarios, the new approach allows millions of scenarios. It also speeds up the computation. \n\n2. The examples shown in the draft are simple but very clear: The 2D and 3D examples shown in the paper are very classical DP problems and the paper make it very clear about how to take advantage of GPU to solve the issues.\n\n3. The experiment results are solid and sound. The results clearly illustrates the advantages of solving the problem using GPU vs CPU.\n\n4. The paper also provides open-source CUDA implementation for reproducibility.", "weaknesses": "1. Limitation of the applicable cases: the authors in the conclusion show that the method is only applicable to some specific type of issues. \nFurther explain in the Questions sections.\n\n2. Some of the notations are not very clear. Further explain in the Questions sections.\n\n3. Minor typos and grammar issues. Further explain in the Questions sections.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)