Pith. sign in

REVIEW 3 major objections 3 minor 30 references

Sample Complexity of Data-driven Multistage Stochastic Programming under Markovian Uncertainty

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

Pith's one-line read The paper claims a polynomial sample bound for multistage stochastic programs under Markovian uncertainty, replacing SAA's exponential dependence on the horizon.

desk verdict Clever and important idea, but the polynomial-in-T sample-complexity claim is not proven — the Bernstein step in Proposition 1 collapses under conditioning. read the letter →

arxiv 2412.19299 v1 pith:CZ2S5PPH submitted 2024-12-26 math.OC

classification math.OC MSC 90C1590C39
keywords multistagestochasticprogrammingsamplecomplexityMarkovianuncertaintyMarkovrecombiningscenariotreeaverageapproximationkernelestimatorcurseofdimensionality
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

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

The reading

This paper proposes the Markov Recombining Scenario Tree (MRST), an approximation to a multistage stochastic program whose random data follow a continuous-state Markov process with an unknown transition distribution. The method needs only two independent historical trajectories: it alternates between them across stages and estimates each conditional expectation with a kernel-weighted average over nearby observed states. The central claim is a non-asymptotic suboptimality bound: with $N \geq \widetilde{O}(T^{p+3}\epsilon^{-p-2})$ samples, the first-stage cost of the MRST solution is within $\epsilon$ of the true optimal value with high probability. Because this bound is polynomial in the horizon $T$, the method is claimed to break the curse of dimensionality that makes sample average approximation need $\widetilde{O}(T^{T}\epsilon^{-2T})$ samples. A sympathetic reader should care because this is a provable polynomial-in-$T$ guarantee for data-driven multistage decisions under Markovian noise.

What carries the argument

The load-bearing object is the Markov Recombining Scenario Tree: a tree assembled from two independent trajectories $\gamma^{(a)}_{[1:N+1]}$ and $\gamma^{(b)}_{[1:N+1]}$, using nodes from the first trajectory at odd stages and from the second at even stages. Conditional expectations are replaced by kernel estimators $\widehat{\mathbb{E}}_{\xi_t}[\cdot] = \sum_i Q_{t+1}(x_t,\gamma^{(a)}_{i+1}) K_h(\xi_t-\gamma^{(a)}_i)/\sum_j K_h(\xi_t-\gamma^{(a)}_j)$ (and the analogous formula with $\gamma^{(b)}$), with a bounded-support kernel $K_h$. The proof's key move is to represent the Markov process as $\tilde{\xi}_{t+1}=\phi(\tilde{\xi}_t,\tilde{\zeta}_{t+1})$ with i.i.d. uniforms $\tilde{\zeta}_t$, so that, conditional on the observed trajectory, the sum inside the estimator is a sum of independent random variables and Bernstein's inequality applies. This converts the non-i.i.d. Markov dependence into a tractable concentration problem, and the recursion in the probability of an empty neighborhood uses the uniform lower bound $f$ on the transition density.

What would settle it

Take a Markov process that satisfies all assumptions except the uniform lower bound on the transition density, for example a transition density that vanishes as the state approaches the boundary of the bounded support, and measure the probability that the kernel estimator finds no neighbor within bandwidth $h$ as $N$ grows; if that probability does not decay exponentially in $N$, the claimed sample complexity fails. Equivalently, simulate Proposition 1 with a truncated Gaussian transition and check directly whether the bias term exceeds $M_{t+1}h+\epsilon$ at the documented rate.

Watch

Extended reading notes

Core claim

On its own terms, the paper establishes Theorem 1: for any fixed $\epsilon \in [0,1]$, the optimal first-stage decision of the MRST problem satisfies $c_1(x_1^*,\xi_1)+\mathbb{E}_{\xi_1}[Q_2(x_1^*,\tilde{\xi}_2)] \leq c_1(\hat{x}_1,\xi_1)+\mathbb{E}_{\xi_1}[Q_2(\hat{x}_1,\tilde{\xi}_2)]+\epsilon$ with high probability once $N \geq \widetilde{O}(T^{p+3}\epsilon^{-p-2})$. The argument proceeds by bounding the error of the kernel estimator for a single conditional expectation (Proposition 1), extending the bound uniformly over the feasible set with an $\eta$-net (Proposition 2), and then propagating the errors backward through the dynamic program (Theorem 2). The final suboptimality bound scales like $\widetilde{O}(T^{(p+3)/(p+2)} N^{-1/(p+2)})$, which is what yields the polynomial sample complexity in $T$.

Load-bearing premise

The transition density must stay bounded below by a positive constant everywhere on a bounded support, so every state has a nontrivial chance of producing a neighbor within the kernel bandwidth; if the density has thin regions or unbounded support, the kernel estimate can return zero and the error bound collapses.

Editorial extensions

If this is right

  • For fixed accuracy $\epsilon$, the sample count grows polynomially with horizon $T$ (power $p+3$), so the method remains tractable for long planning horizons where SAA does not.
  • The MRST value functions are defined on observed states, and the same kernel estimator extends the policy to unseen data, so the output is an implementable decision rule rather than a tree-dependent policy.
  • The suboptimality bound degrades only polynomially in the state dimension $p$, and the bandwidth formula $h^\star$ derived from the bound gives an explicit tuning recipe.
  • Under the paper's regularity assumptions, two trajectories suffice for a non-asymptotic guarantee, so the method applies where history data is scarce.

Reading between the lines

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

  • If the bound holds, MRST-style recombining trees could become a default fallback for data-limited multistage optimization, for instance by feeding MRST stage problems into SDDP-type solvers.
  • The uniform lower bound on the transition density is probably not satisfied by the Gaussian noise used in the LQG experiments unless truncated; testing on a state-dependent or unbounded transition law would reveal how far the guarantee extends beyond Assumption A3.
  • A natural next question is whether the polynomial-in-$T$ rate can be improved by using more than two trajectories or by replacing the kernel step with a regression estimator, since the $p$-dependence looks like a smoothing cost rather than an inherent horizon curse.
  • The alternating-trajectory construction implicitly exploits stationarity of the Markov chain; re-proving the bound under mixing coefficients instead of a uniform density lower bound would connect these guarantees to a wider class of dependent data.
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 / 3 minor

Summary. The paper proposes a Markov Recombining Scenario Tree (MRST) method for multistage stochastic programming under an unknown continuous-state Markov process. The method uses two independent historical trajectories and kernel-weighted conditional expectations to build an approximate dynamic program. The main claimed contribution is a non-asymptotic sample complexity of eO(T^{p+3} eps^{-p-2}) for an eps-optimal first-stage solution, in contrast to the exponential-in-T complexity of SAA. The theoretical argument is built on Proposition 1, a concentration inequality for a kernel estimator along a single Markovian trajectory; this proposition is then used to derive uniform generalization bounds, a recursive generalization bound for the approximate value functions, and a suboptimality bound. Numerical experiments on a Linear Quadratic Gaussian control problem compare MRST with SAA.

Significance. If the main theorem were correct, the result would be significant: it would provide the first non-asymptotic sample-complexity guarantee for data-driven multistage stochastic programming under Markovian uncertainty that is polynomial in the horizon, and the proposed estimator is simple and computationally attractive. The paper is well organized, the regularity assumptions are explicit, and the derivations are not circular: the bandwidth h and net width eta are proof parameters rather than quantities fitted to data. The claimed eO(T^{p+3} eps^{-p-2}) bound would address a well-known open challenge. However, the central concentration proof contains a load-bearing error, and the numerical experiments are run under assumptions that the theoretical part does not cover. As it stands, the main claim is not established.

major comments (3)
  1. [Section 3, Proposition 1, Eqs. (21)-(23) and (35)] The proof of Proposition 1 is invalid at the Bernstein step. The summands in Eq. (25) are D_i = E_{gamma_i}[Q(x_t, phi(gamma_i, zeta_{i+1}))] - Q(x_t, phi(gamma_i, zeta_{i+1})) with kernel weights. The proof conditions on gamma_i = gamma_i for all i in [N] and then claims that the D_i are independent because zeta_2, ..., zeta_{N+1} are i.i.d. But conditioning on the whole trajectory gamma_1, ..., gamma_N fixes gamma_2, ..., gamma_N, and hence, through gamma_{i+1} = phi(gamma_i, zeta_{i+1}), it fixes zeta_{i+1} for i < N; only zeta_{N+1} remains random. The conditioned expression is therefore an affine function of one random variable plus the deterministic offset sum_{i<N} w_i (E_{gamma_i}[Q] - Q(gamma_i, gamma_{i+1})), not a sum of N independent centered random variables. Consequently, Bernstein's inequality cannot be applied as in Eqs. (27)-(29). Since Corollary 1, Proposition 2, Theorem 2, Corollary 2, and Theorem 1 all invoke Proposition 1, the polynomial sample complexity claim is not supported by the submitted proof.
  2. [Section 3, Assumption A3, Eqs. (21)-(23), (35)] The proof uses the estimate integral_{Xi_h(gamma)} f(xi | gamma') dxi >= C_p h^p f for every gamma' in Xi, where Xi_h(gamma) = Xi cap B(gamma, h). A uniform lower bound on the transition density f(xi | gamma) over the bounded set Xi does not imply this estimate: for states gamma near the boundary of Xi, the intersection Xi cap B(gamma, h) can have volume much smaller than C_p h^p. The same issue affects the exponential moment bound in Eq. (35). The paper would need an additional assumption, for example a uniform lower bound on the transition mass into every h-ball or a restriction that gamma is bounded away from the boundary of Xi. As written, Proposition 1 fails even setting aside the independence problem in Eqs. (25)-(29).
  3. [Section 4, Eqs. (55)-(56), Figure 2] The numerical experiments do not operate under the assumptions of Theorem 1. The LQG specification uses additive Gaussian noise, so the state process has unbounded support and violates the bounded-support and uniform lower-density conditions in Assumption A3, and the quadratic cost is not scaled to lie in [0,1] as required by Assumption A2. No truncation or rescaling is reported. The favorable numerical comparison with SAA therefore cannot be used as evidence for the theoretical sample-complexity claim, and the experiments do not test the theorem's hypotheses.
minor comments (3)
  1. [Theorem 1] Theorem 1 states 'with high probability' without specifying the failure probability; it should state explicitly that for every delta in (0,1], the claimed bound holds with probability at least 1-delta when N is at least the stated order.
  2. [Assumption A3] The same symbol f is used for the transition density and for its lower bound in the statement '0 < f <= f(xi | gamma)'. This is confusing and should be disambiguated, for example by writing f_lower for the lower bound.
  3. [Figure 2 and Section 4] The x-axis of Figure 2 is labeled as the number of nodes in the scenario tree, but MRST and SAA consume data in different ways: MRST uses only two trajectories of length N+1, while the SAA tree size is not the same as a data sample count. The comparison would be clearer if the paper reported the actual number of historical observations used by each method.

Circularity Check

0 steps flagged · score 1.0 of 10

No significant circularity: the MRST sample-complexity bound is derived from explicit kernel-estimator generalization bounds with analytically optimized bandwidth, not from fitted inputs or self-citation.

full rationale

The derivation chain is self-contained. Proposition 1 splits the error of the kernel estimator (5)-(6) into a bias term (a) bounded by Mt+1h plus a no-neighbor event (Eqs. (18)-(24)) and a deviation term (b) bounded by Bernstein's inequality (Eqs. (25)-(29)); Corollary 1 inverts the exponential tail; Proposition 2 adds an eta-net union bound; Theorem 2 accumulates the bounds backward through the dynamic program; Corollary 2 derives the suboptimality bound and then optimizes h analytically (Eqs. (49)-(54)) to obtain eO(T^{p+3} epsilon^{-p-2}). None of the bounds is obtained by fitting a parameter to the data being predicted: h and eta are proof parameters chosen by minimization, and the sample complexity is the solution of the resulting inequality for N. The paper's only self-citation, [14] (Park, Jia, Hanasusanto), appears in Remark 3 as one entry in a list of SDDP variants and is explicitly labeled 'beyond the scope of this paper', so it is not load-bearing. The external references used in the proof (Kallenberg's Markov representation, Vershynin's Bernstein inequality, and Shapiro-Nemirovski's net cardinality) are standard mathematical facts, not author-specific results that smuggle in the conclusion. A possible gap in Proposition 1's conditioning argument, whereby conditioning on e_gamma_i = gamma_i for all i makes all but the last summand deterministic so the stated Bernstein application may be unjustified, is a correctness or rigor concern, not circularity, because accepting it would not make Theorem 1 equivalent to an input or to a fitted value.

Assumptions & free parameters 2 free parameters · 7 assumptions · 0 invented entities

The central claim depends on the standard regularity assumptions A1-A5, on the Markov representation, and on the eta-net cardinality bound. The tunable parameters are the kernel bandwidth h and net width eta. No new physical or mathematical entities are introduced; the MRST is an algorithmic construction rather than an invented entity.

free parameters (2)
  • bandwidth h = h* proportional to (T/N)^{1/(p+2)} in Eq. (53)-(54); not specified in experiments
    The kernel estimator in (5)-(6) depends on the bandwidth h. The theoretical bound is optimized over h, but the numerical experiments do not state how h was selected, leaving a potential tuning parameter.
  • net width eta = not specified
    The uniform generalization bound in Proposition 2 uses an eta-net of the feasible decision regions. The final bound depends on eta through an L*eta term and logarithmic factors, but no specific value is recommended or reported in experiments.
assumptions (7)
  • domain assumption A1: feasible regions X_t have finite diameter D_t uniformly over states
    Used to bound the size of eta-nets in Proposition 2 and Theorem 2; a standard regularity condition in MSP complexity analysis.
  • domain assumption A2: value functions Q_t are bounded in [0,1] and Lipschitz in the decision vector
    Used throughout the proof for Lipschitz constants L_t and to bound value functions and their variances. The paper argues scaling can enforce the [0,1] range, but Lipschitzness of the true value function is a nontrivial property of the problem data.
  • domain assumption A3: time-homogeneous Markov process with bounded support and transition density lower bounded by f>0
    The lower density bound is used in Proposition 1 to guarantee that a ball of radius h has probability at least C_p h^p f of containing each trajectory point. This is a strong assumption and is violated by unbounded supports such as Gaussian noise.
  • domain assumption A4: conditional expectations E_gamma[Q] and E_gamma[Q^2] are Lipschitz in the conditioning state with constant M_t
    Used to bound the bias term M_t h and the variance Lipschitz constant in Lemma 3. It requires the future value function's conditional mean to vary smoothly with the current state, which is not guaranteed for arbitrary costs and dynamics.
  • domain assumption A5: kernel function k is non-increasing and bounded above
    Standard for Nadaraya-Watson type estimators; used to bound weights in Proposition 1.
  • standard math Definition 1: every Markov process has a representation xi_{t+1} = phi(xi_t, zeta_{t+1}) with i.i.d. uniforms zeta
    Kallenberg's Markov representation is invoked in the proof of Proposition 1 to decouple the innovations. The proof's conditioning step is where this representation is misused.
  • standard math Cardinality of eta-net of a bounded set is O(1)(D/eta)^d
    Taken from Shapiro and Nemirovski and used in Proposition 2 and Theorem 2 to apply union bounds over discretized decision sets.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Sample Complexity of Data-driven Multistage Stochastic Programming under Markovian Uncertainty." pith.science (2026). https://pith.science/paper/CZ2S5PPH

@misc{pith2026241219299,
  author       = {Pith},
  title        = {Pith review of: Sample Complexity of Data-driven Multistage Stochastic Programming under Markovian Uncertainty},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/CZ2S5PPH}},
  note         = {Machine review of arXiv:2412.19299}
}
abstract

This work is motivated by the challenges of applying the sample average approximation (SAA) method to multistage stochastic programming with an unknown continuous-state Markov process. While SAA is widely used in static and two-stage stochastic optimization, it becomes computationally intractable in general multistage settings as the time horizon $T$ increases. Indeed, the number of samples required to obtain a reasonably accurate solution grows exponentially$\text{ -- }$a phenomenon known as the curse of dimensionality with respect to the time horizon. To overcome this limitation, we propose a novel data-driven approach, the Markov Recombining Scenario Tree (MRST) method, which constructs an approximate problem using only two independent trajectories of historical data. Our analysis demonstrates that the MRST method achieves polynomial sample complexity in $T$, providing a more efficient alternative to SAA. Numerical experiments on the Linear Quadratic Gaussian problem show that MRST outperforms SAA, addressing the curse of dimensionality.

Figures

Figures reproduced from arXiv: 2412.19299 by the authors.

Figure 1
Figure 1. Scenario tree using two historical trajectories when T = 5. these in-sample policies are not directly implementable in the true environment, where the data process may differ from the in-sample data. Nevertheless, our approach allows us to derive an implementable policy xbt(ξ[1:t] ) for any unseen data. Given the history ξ[1:t] and the approximate value function Qbt+1(·), we can compute the out-of-sample policy xbt(… view at source ↗
Figure 2
Figure 2. Comparison of optimality gaps between SAA and MRST: solid lines represent the out-of-sample suboptimality for MRST (blue) and SAA (yellow), while dashed lines represent the in-sample suboptimality. 26 [PITH_FULL_IMAGE:figures/full_fig_p026_2.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

30 extracted references · 27 canonical work pages

  1. [1]

    J., and Parija, G

    Ahmed, S., King, A. J., and Parija, G. (2003). A multi-stage stochastic integer programming approach for capacity expansion under uncertainty. Journal of Global Optimization , 26:3–24

  2. [2]

    Bellman, R. (1957). Dynamic Programming. Princeton University Press

  3. [3]

    F., Cen, Z., and Christel, T

    Bonnans, J. F., Cen, Z., and Christel, T. (2012). Energy contracts management by stochastic programming techniques. Annals of Operations Research, 200:199–222

  4. [4]

    Markov Chain-based Policies for Multi-stage Stochastic Integer Linear Programming with an Application to Disaster Relief Logistics

    Castro, M. P., Bodur, M., and Song, Y. (2022). Markov chain-based policies for multi-stage stochastic integer linear programming with an application to disaster relief logistics. arXiv preprint arXiv:2207.14779

  5. [5]

    and Morton, D

    Duque, D. and Morton, D. P. (2020). Distributionally robust stochastic dual dynamic program- ming. SIAM Journal on Optimization , 30(4):2841–2865. 29

  6. [6]

    computational complexity of stochastic programming problems

    Hanasusanto, G. A., Kuhn, D., and Wiesemann, W. (2016). A comment on “computational complexity of stochastic programming problems”. Mathematical Programming, 159:557–569

  7. [7]

    Huang, J., Zhou, K., and Guan, Y. (2017). A study of distributionally robust multistage stochastic optimization. arXiv preprint arXiv:1708.07930

  8. [8]

    and Morton, D

    Infanger, G. and Morton, D. P. (1996). Cut sharing for multistage stochastic linear programs with interstage dependency. Mathematical Programming, 75(2):241–256

Show all 30 references
  1. [9]

    and Li, S

    Jiang, J. and Li, S. (2021). On complexity of multistage stochastic programs under heavy tailed distributions. Operations Research Letters, 49(2):265–269

  2. [10]

    Kallenberg, O. (1997). Foundations of modern probability, volume 2. Springer

  3. [11]

    J., Shapiro, A., and Homem-de Mello, T

    Kleywegt, A. J., Shapiro, A., and Homem-de Mello, T. (2002). The sample average approxima- tion method for stochastic discrete optimization. SIAM Journal on optimization , 12(2):479–502

  4. [12]

    and Shapiro, A

    L¨ ohndorf, N. and Shapiro, A. (2019). Modeling time-dependent randomness in stochastic dual dynamic programming. European Journal of Operational Research, 273(2):650–661

  5. [13]

    Lozier, D. W. (2003). Nist digital library of mathematical functions. Annals of Mathematics and Artificial Intelligence , 38:105–119

  6. [14]

    Park, H., Jia, Z., and Hanasusanto, G. A. (2022). Data-driven stochastic dual dynamic pro- gramming: Performance guarantees and regularization schemes. Available at Optimization On- line

  7. [15]

    Pereira, M. V. and Pinto, L. M. (1991). Multi-stage stochastic optimization applied to energy planning. Mathematical Programming, 52(1):359–375

  8. [16]

    Philpott, A. B. and De Matos, V. L. (2012). Dynamic sampling algorithms for multi-stage stochastic programs with risk aversion.European Journal of operational research, 218(2):470–483

  9. [17]

    B., de Matos, V

    Philpott, A. B., de Matos, V. L., and Kapelevich, L. (2018). Distributionally robust sddp. Computational Management Science , 15:431–454

  10. [18]

    Reaiche, M. (2016). A note on sample complexity of multistage stochastic programs.Operations Research Letters, 44(4):430–435. 30

  11. [19]

    Robbins, H. (1955). A remark on stirling’s formula. The American mathematical monthly , 62(1):26–29

  12. [20]

    and Petersen, I

    Shaiju, A. and Petersen, I. R. (2008). Formulas for discrete time lqr, lqg, leqg and minimax lqg optimal control problems. IF AC Proceedings Volumes, 41(2):8773–8778

  13. [21]

    Shapiro, A. (2006). On complexity of multistage stochastic programs. Operations Research Letters, 34(1):1–8

  14. [22]

    Shapiro, A. (2011). Analysis of stochastic dual dynamic programming method. European Journal of Operational Research, 209(1):63–72

  15. [23]

    Shapiro, A., Dentcheva, D., and Ruszczynski, A. (2021). Lectures on stochastic programming: modeling and theory . SIAM

  16. [24]

    and Nemirovski, A

    Shapiro, A. and Nemirovski, A. (2005). On complexity of stochastic programming problems. Continuous optimization: Current trends and modern applications , pages 111–146

  17. [25]

    P., and Soares, M

    Shapiro, A., Tekaya, W., da Costa, J. P., and Soares, M. P. (2013). Risk neutral and risk averse stochastic dual dynamic programming method. European journal of operational research, 224(2):375–391

  18. [26]

    Silva, T., Vallad˜ ao, D., and Homem-de Mello, T. (2021). A data-driven approach for a class of stochastic dynamic optimization problems. Computational Optimization and Applications , 80:687–729

  19. [27]

    and Wasserman, L

    Tibshirani, R. and Wasserman, L. (2013). Nonparametric regression. Statistical Machine Learning, Spring

  20. [28]

    Vershynin, R. (2018). High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press

  21. [29]

    and Sun, X

    Zhang, S. and Sun, X. A. (2022). Stochastic dual dynamic programming for multistage stochas- tic mixed-integer nonlinear optimization. Mathematical Programming, 196(1):935–985

  22. [30]

    Zou, J., Ahmed, S., and Sun, X. A. (2019). Stochastic dual dynamic integer programming. Mathematical Programming, 175:461–502. 31

Pith tools

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