Polynomial-Time Algorithm for Optimal Stopping with Fixed Accuracy

1. Introduction

We consider the problem of discrete-time optimal stopping (OS): a decision maker observes a stochastically evolving information process $\mathbf{Y}=\{Y_1,\dots ,Y_T\}$, where we assume that (for all t) $Y_t$ is a D-dimensional vector for some $D\ge 1$. Let $\mathcal{F}=\{{\mathcal{F}}_t,t=1,\dots ,T\}$ be the natural filtration generated by $\mathbf{Y}$, that is, ${\mathcal{F}}_t\triangleq \sigma (Y_1,\dots ,Y_t)$. Let $\{g_t:{\mathbb{R}}^{D\times t}\to \mathbb{R},t=1,\dots ,T\}$ be a sequence of measurable reward functions mapping trajectories $Y_{[t]}$ to reals, with $Y_{[t]}=(Y_1,\dots ,Y_t)$. For notational simplicity, let $Z_t\triangleq g_t(Y_{[t]}),$ and denote the corresponding stochastic process by $\mathbf{Z}$. Throughout this paper, we assume $Z_t$ is integrable for all t, and the problem of interest is to compute $\text{OPT}={\sup }_{\tau \in \mathcal{T}}E[Z_{\tau }]$, where $\mathcal{T}$ denotes the set of all $\mathcal{F}$-adapted stopping times.

Problems of this kind find an array of applications across areas including operations research and management science, finance, economics, statistics, and computer science. In practice, the problems are often high-dimensional, with long horizons and non-Markovian/path-dependent dynamics, posing significant computational challenges. These challenges have generated substantial research across multiple academic communities and led to a vast literature on approximation algorithms and heuristics. Here, we briefly describe some of the main past approaches to OS, and defer a more detailed discussion to our literature review in Section 1.2.

The canonical approach to OS is to formulate it as a stochastic dynamic program (DP). The associated Bellman equation takes the form $V_t=\max (Z_t,E[V_{t+1}|{\mathcal{F}}_t]),t=1,\dots ,T-1$, and $V_T=Z_T.$ Here, $V_t$ is the reward-to-go function in period t, and we note that $V_t$ is ${\mathcal{F}}_t$-measurable for all t. Then $\text{OPT}=E[V_1].$ Naively implementing this recursion, however, takes time exponential in the dimension and/or time horizon, a computational challenge known as the curse of dimensionality. Since the late 1990s, a major line of research has sought to overcome this challenge through approximate dynamic programming (ADP) methods (Tsitsiklis and Van Roy 1999, Longstaff and Schwartz 2001). Here one uses a judicious choice of basis functions and/or random sampling to approximate the Bellman equation and yield tractable algorithms, although typically having an optimality gap that depends on the quality of these basis functions (and that cannot be made arbitrarily small without incurring an exponential complexity or making extra assumptions).

Another popular approach builds on a dual formulation for OS, pioneered in the seminal work of Davis and Karatzas (1994). They establish a strong duality: $\text{OPT}={\inf }_{\mathbf{M}\in {\mathcal{M}}^0}E[{\max }_{t=1,\dots ,T}(Z_t-M_t)]$, where ${\mathcal{M}}^0$ denotes the set of mean-zero martingales adapted to $\mathcal{F}$. Moreover, they showed that the Doob–Meyer decomposition of the value functions yields an optimal martingale ${\mathbf{M}}^{*}$. Although not offering a resolution to the curse of dimensionality (as computing ${\mathbf{M}}^{*}$ requires prior knowledge of the value functions themselves), this duality theory led to substantial algorithmic progress that complements ADP methods.

In this work, we provide new theoretical insights into OS by introducing a novel recursive expansion representation of the OS value, which leads to a polynomial-time algorithm with provable approximation guarantees. Below we summarize our main contributions.

1.1. Main Contributions

1.1.1. Novel Expansion Theory.

The standard DP formulation essentially relies on backward induction (in time) to set up the recursive representation of opt. We develop a novel expansion representation of opt that builds on fundamentally different logic. More precisely, the first term of our expansion coincides with the expected pathwise maximum $E[{\max }_{t=1,\dots T}{Z}_t]$, that is, the hindsight optimal value. The second term is an expectation that can be interpreted as a type of regret, and subsequent terms are defined recursively as the expectation of a maximum that can be interpreted as a higher-order notion of regret. This expansion representation has an intimate connection to the dual martingale formulation and yields an optimal dual martingale different from ${\mathbf{M}}^{*}.$ Under certain natural regularity conditions on the reward process (e.g., Assumption 1 on bounded rewards, or Assumption A.1 in the appendix regarding the bounded squared coefficient of variation of the maximum reward), we show that truncating our expansion after $O\left(\frac{1}{ϵ}\right)$ terms yields an $ϵ$-approximation to opt, with those terms of the expansion requiring a level of nesting (of conditional expectations) scaling only as $O\left(\frac{1}{ϵ}\right)$ (independent of the time horizon or dimension). We believe this representation complements the existing literature and offers a novel perspective on the theory of OS.

1.1.2. Polynomial-Time Algorithm.

We formalize an algorithm that implements (i.e., approximately computes) an appropriate truncation of our expansion representation via Monte Carlo simulation, given access to an efficient simulator of $\mathbf{Y}$. Using this algorithm, we implement an $ϵ$-optimal policy in time $T^{O\left(\frac{1}{ϵ}\right)}$, for any fixed $ϵ\in (0,1)$. More precisely, the optimal dual martingale associated with our approach has an analogous expansion representation. Our stopping policy truncates this expansion at depth $O\left(\frac{1}{ϵ}\right)$, uses Monte Carlo simulation to approximate the value of this truncated martingale at each time period, and stops the first time the reward process is sufficiently close to this approximated martingale. As simulating the kth term of the expansion requires implementing depth k nested conditional expectations, this requires a simulation and computational complexity scaling as $T^{O\left(\frac{1}{ϵ}\right)}$. Because many other approaches (including DP and naive fully nested simulation) have a worst-case complexity growing exponentially with the horizon and/or dimension, our results demonstrate an exponential improvement and serve as a proof of concept that general high-dimensional path-dependent OS problems are polynomially approximable. To be clear, here polynomial time is meant in the fixed-accuracy sense: for each fixed $ϵ$, the exponent of T is a constant (depending on $ϵ$). Our algorithms have an exponential (not polynomial) dependence on $\frac{1}{ϵ}$ itself.

1.2. Related Literature

There is a vast literature on computational OS, and one of the most popular approaches is regression/ADP. Here one fixes a family of basis functions to approximate the DP value functions. Seminal papers in this area include Tsitsiklis and Van Roy (2001) and Longstaff and Schwartz (2001). There was much subsequent work on, for example, nonparametric regression, smoothing spline regression, recursive kernel regression, integer programming, reinforcement learning, multilevel Monte Carlo, and (deep) neural networks (see Lai and Wong 2004; Egloff 2005; Kohler 2008; Belomestny et al. 2015, 2020; Bayer et al. 2021; Becker et al. 2021; Ciocan and Mišić 2022; Sturt 2023; Zhou et al. 2023). These methods mostly focus on achieving empirical success rather than ensuring theoretical guarantees, with their performance relying on how well the choice of basis functions can approximate the true value functions. Those works with theoretical guarantees often require additional continuity assumptions, and/or scale exponentially in the time horizon. Most existing theoretical analyses for these methods focus on their convergence to the best approximation within the fixed set of basis functions (Clément et al. 2002, Glasserman and Yu 2004, Stentoft 2004, Belomestny 2011, Bouchard and Warin 2012, Bezerra et al. 2020) and do not have guarantees comparable to our own.

Building on the seminal work of Davis and Karatzas (1994), significant progress was made simultaneously by Haugh and Kogan (2004) and Rogers (2002) in their formulation of a dual methodology for OS. Instead of finding an OS time, it formulates and solves a dual optimization problem over the set of adapted martingales. Other dual representations were subsequently discovered (e.g., Jamshidian 2007), and the methodology was extended to more general control problems (Brown et al. 2010). We note that closely related dual formulations and methodologies have also appeared in the study of multistage stochastic optimization (Rockafellar and Wets 1991). Simulation approaches via dual formulation have since led to substantial algorithmic progress (see Andersen and Broadie 2004; Kolodko and Schoenmakers 2004, 2006; Belomestny and Milstein 2006; Chen and Glasserman 2007; Desai et al. 2012; Belomestny et al. 2013, 2019; Christensen 2014; Lelong 2018; Ibáñez and Velasco 2020). To achieve computational efficiency, many of these algorithms adopt ideas similar to those of ADP approaches, for example, fixing a family of basis functions/martingales to approximate the optimal dual martingale. Consequently, the quality of the bounds derived using these approaches typically depends on the expressiveness of the set of basis functions/martingales, or the accuracy of the initial estimation of the value functions. In addition, those works that do provide rigorous approximation guarantees typically require a level of nesting (of conditional expectations) scaling with the time horizon (see Kolodko and Schoenmakers 2006, Chen and Glasserman 2007).

More recently, new machine learning techniques, such as decision trees and deep neural networks, and ideas from robust optimization were introduced to solve the problem of high-dimensional OS and American option pricing (see Kohler et al. 2010, Becker et al. 2019, Becker et al. 2021, Ciocan and Mišić 2022, Sturt 2023). These approaches exhibit excellent performance on a number of benchmarks, but generally do not have strong theoretical guarantees independent of the quality of the basis functions/free from additional continuity assumptions.

1.3. Organization

The rest of this paper is organized as follows. In Section 2, we describe our expansion representation, establish its connection to the martingale duality theory, and provide associated error bounds. In Section 3, we introduce our expansion-inspired algorithm and state and prove our main algorithmic results. We present some closing thoughts and directions for future research in Section 4. We include an appendix, which contains several auxiliary technical lemmas used in our proof, along with a generalization to unbounded rewards. We also include an Online Companion, in which we prove that the well-studied Bermudan max call falls within our theoretical framework and present the results of some numerical experiments (mostly to illustrate the practical limitations of our approach).

2. Novel Expansion Representation for opt

In this section, we establish our expansion representation of opt, draw connections to the martingale duality theory, and analyze its convergence rate.

2.1. Novel Expansion Representation for opt

We begin by giving the simple intuition behind the expansion. We wish to compute $\text{OPT}={\sup }_{\tau \in \mathcal{T}}E[Z_{\tau }].$ Trivially, ${\sup }_{\tau \in \mathcal{T}}E[Z_{\tau }]\le E[{\max }_{t=1,\dots ,T}{Z}_t].$ We next turn this straightforward bound into an equality by compensating with a remainder term, which is characterized by a new OS problem. We introduce a specific martingale $\{E[{\max }_{i=1,\dots ,T}{Z}_i|{\mathcal{F}}_t],\hspace{0.17em}t=1,\dots ,T\},$ that is, the Doob martingale of ${\max }_{t=1,\dots ,T}{Z}_t,$ adapted to filtration $\mathcal{F}.$ By the optional stopping theorem, for any $\tau \in \mathcal{T}$, it holds that $E[E[{\max }_{t=1,\dots ,T}{Z}_t|{\mathcal{F}}_{\tau }]]=E[{\max }_{t=1,\dots ,T}{Z}_t],$ and thus

Taking supremum on both sides, we conclude that

For $t=1,\dots ,T$, let $Z_t^1\triangleq Z_t,Z_t^2\triangleq Z_t^1-E[{\max }_{i=1,\dots ,T}{Z}_i^1|{\mathcal{F}}_t]$ (with probability one (w.p.1) nonpositive). Then (1) becomes

There are many bounds for the gap $E[{\max }_{t=1,\dots ,T}{Z}_t^1]-\text{OPT}$ in the literature on so-called prophet inequalities (see Hill and Kertz 1983, 1992). The novelty of our approach will be to “recurse” this logic on the problem ${\sup }_{\tau \in \mathcal{T}}E[Z_{\tau }^2]$ itself, and subsequently a sequence of such problems defined recursively. More formally, let $Z_t^{k+1}\triangleq Z_t^k-E[{\max }_{i=1,\dots ,T}{Z}_i^k|{\mathcal{F}}_t]$ for all $k\ge 1$ and $t\in \{1,\dots ,T\}$ (denoting the corresponding stochastic process by ${\mathbf{Z}}^k$).

It follows from the basic properties of conditional expectation and a straightforward induction (the details of which we omit) that $Z_t^k\le 0$ for all $t\in \{1,\dots ,T\}$ and $k\ge 2$. Then our first main result is as follows.

Before proving Theorem 1, we first state and prove some auxiliary results for ${\mathbf{Z}}^k$.

With Lemma 1 in hand, we would be done with the proof of Theorem 1 if we could show that ${\lim }_{k\to \infty }{\sup }_{\tau \in \mathcal{T}}E[Z_{\tau }^{k+1}]=0.$ We now prove that this is indeed the case.

Let us more explicitly give the first few terms. For $k\ge 1$, let $L_k\triangleq E[{\max }_{t=1,\dots ,T}{Z}_t^k]$. Then

Note that the first term, $L_1$, is the only positive term in the expansion, which corresponds to the obvious upper bound. The term $L_2$ has a “regret” interpretation, as it is the expectation of the maximum of T terms, with term t intuitively the reward for stopping at time t minus your best guess of the all-time maximum given the information up to time t (i.e., the regret of stopping at time t). Later terms are the expectations of similar explicit suprema with analogous interpretations in terms of certain notions of (negative) “higher-order” regret, each of which can be computed by simulation.

2.2. Connection to Martingale Duality

As a byproduct of the expansion representation, we derive a novel optimal dual martingale. Let $S\triangleq {\sum }_{k=1}^{\infty }{\max }_{t=1,\dots ,T}{Z}_t^k$. Integrability of S is easily seen to follow from Theorem 1 and monotone convergence (as opt is finite). Let MAR denote the Doob martingale such that (s.t.) ${\mathit{\text{MAR}}}_t=E[S|{\mathcal{F}}_t]$, $t=1,\dots ,T$. We now prove that $\mathbf{M}\mathbf{A}{\mathbf{R}}^0\triangleq \mathbf{M}\mathbf{A}\mathbf{R}-E[S]$ is an optimal dual martingale for OS, and begin by proving the following result.

Because $P({\max }_{t=1,\dots ,T}(Z_t-{\mathit{\text{MAR}}}_t)=0)=1$ implies $E[{\max }_{t=1,\dots ,T}(Z_t-({\mathit{\text{MAR}}}_t-\text{OPT}))]=\text{OPT}$, Lemma 2 immediately implies that $\mathbf{M}\mathbf{A}{\mathbf{R}}^0$ is indeed an optimal dual martingale in the sense of Davis and Karatzas (1994).

We now discuss a more refined notion of optimality to delineate a key difference between $\mathbf{M}\mathbf{A}{\mathbf{R}}^0$ and the standard optimal dual martingale ${\mathbf{M}}^{*}$ arising from the Doob–Meyer decomposition of the value functions. First, we recall the notion of i-sure optimality, as formalized in Schoenmakers et al. (2013) (our definition also adjusts for the different time indexing used there).

By Lemma 2, ${\max }_{j=1,\dots ,T}(Z_j-MA{R}_j^0+MA{R}_1^0)=MA{R}_1^0$. Thus, $\mathbf{M}\mathbf{A}{\mathbf{R}}^0$ would be 0-surely optimal if w.p.1 $MA{R}_1^0=V_1$ (the first reward-to-go function). As by definition $MA{R}_1^0=E[{\sum }_{k=1}^{\infty }{\max }_{t=1,\dots ,T}{Z}_t^k|{\mathcal{F}}_1]$, the desired equality follows from a proof essentially identical to our proof of Theorem 1, but with all expectations replaced by conditional expectations taken with respect to (w.r.t.) ${\mathcal{F}}_1$, and we omit the details. We conclude the following.

We note that (as shown in Schoenmakers et al. 2013) 0-sure optimality is a stronger notion than optimality in the sense of Davis and Karatzas (1994), as it requires optimality in a pathwise sense (for the zero-mean martingale $\{M_t-M_1,t=1,\dots ,T\}$), not just in expectation. It is interesting to point out that the standard optimal dual martingale ${\mathbf{M}}^{*}$ of Davis and Karatzas (1994) (arising from the Doob–Meyer decomposition of the value functions) is not only 0-surely optimal, but also surely optimal (see Schoenmakers et al. 2013). We now show that our optimal dual martingale $\mathbf{M}\mathbf{A}{\mathbf{R}}^0$ does not satisfy this stronger sure optimality property, creating a clear delineation between $\mathbf{M}\mathbf{A}{\mathbf{R}}^0$ and ${\mathbf{M}}^{*}$, and defer the proof to the appendix, Section A.1.

Next, we point out that ${\tau }_{\infty }$ can be interpreted as a “dual-threshold” stopping time w.r.t. MAR, and is in fact an optimal stopping time. Recall that ${\tau }_{\infty }=\min \{t\in \{1,\dots ,T\}:Z_t^{\infty }=0\}$, and hence, by (6), ${\tau }_{\infty }=\min \{t\in \{1,\dots ,T\}:Z_t-{\mathit{\text{MAR}}}_t=0\},$ that is, ${\tau }_{\infty }$ can indeed be interpreted as a dual-threshold stopping time.

The polynomial-time stopping policy we define in Section 3 essentially implements the stopping time ${\tau }_{\infty }$ defined in Observation 1, but with $Z^{\infty }$ replaced by the approximation $Z^k$ (for $k=O\left(\frac{1}{ϵ}\right)$) in the definition ${\tau }_{\infty }=\min \{t\in \{1,\dots ,T\}:Z_t^{\infty }=0\}$ (and “equals zero” replaced by “at least $-\frac{1}{k-1}$”), where $Z^k$ is itself approximated by a depth k nested Monte Carlo simulation.

Let us close this section by comparing our approach to other duality-based approaches appearing in the literature. Essentially all of these previous works either (1) assume that one is given a good approximation to the value function (Belomestny et al. 2009) or a set of basis functions whose span contains a good approximation to the value function (Desai et al. 2012) and use these to construct a good candidate martingale; (2) take a primal-dual approach by starting with a candidate (in general suboptimal) stopping time and generating a (in general suboptimal) martingale by considering the martingale part of the Doob–Meyer decomposition of the value process generated by the stopping time (Andersen and Broadie 2004); (3) construct approximately optimal martingales using neural networks (Ye and Wong 2025); or (4) construct a sequence of recursively defined martingales that (after T iterations) converges to the optimal dual martingale (Kolodko and Schoenmakers 2006). Our approach is conceptually most similar to that of Kolodko and Schoenmakers (2006), directly constructing a sequence of recursively defined martingales without needing to be given approximations to the value functions or an associated basis. The key difference between our work and past such recursive constructions is that our sequence of martingales provably converges to optimality at rate $O\left(\frac{1}{k}\right)$, whereas previous such constructions (although exhibiting good empirical performance) had no such theoretical guarantees on the optimality gap of the iterates (beyond achieving optimality after the full T iterations, i.e., any theoretical guarantees required a deep level of nesting for long-horizon problems).

2.3. Approximation Guarantees When Truncating the Expansion

The power of Theorem 1 is that it allows for rigorous approximation guarantees when the infinite sum is truncated, under a natural boundedness assumption on $\mathbf{Z}.$

Let $E_k\triangleq {\sum }_{i=1}^k{L}_i$ denote the finite-sum truncation of our expansion after the first k terms.

Thus, truncating our expansion after k terms yields an absolute error at most $\frac{U}{k}$, independent of the distribution of the underlying stochastic processes $\mathbf{Y}$ and $\mathbf{Z}$ or length of the time horizon, requiring only that the reward functions are bounded in [0, U]. Analogous results also hold when $\mathbf{Z}$ is allowed to be negative (and is uniformly bounded both above and below), although we do not formalize that extension here.

To prove Theorem 2, we state and prove a stronger pathwise guarantee, relating ${\max }_{t=1,\dots ,T}{Z}_t^k$ to ${\max }_{t=1,\dots ,T}{Z}_t^{\infty }$ (which by Lemma 2 and (6) equals zero). By the logic described at the end of Section 2.2, this will ultimately enable us to construct provably accurate and polynomial-time stopping policies by (intuitively) stopping the first time that $Z_t^k\ge -\frac{U}{k-1}$ (for $k=O\left(\frac{U}{ϵ}\right)$), instead of the truly optimal first time that $Z_t^{\infty }\ge 0$.

We next present a lower bound showing that Theorem 2 is tight up to an absolute multiplicative constant factor (asymptotically in k). Namely, there is a precise sense in which our bound on the rate of convergence of the stated expansion cannot be substantially improved. We defer the proof to the appendix, Section A.2.

Thus, under Assumption 1, truncating our expansion after k terms yields an error that is at most $\frac{1}{k}$ and (for large k) at least (asymptotically) $\frac{e^{-1}}{k}$ (in the worst case).

In this section, we established the convergence rate of our expansion under a boundedness assumption on the reward process (Assumption 1). This condition is not essential for our main results, but rather assumed to enhance readability. Essentially the same set of convergence guarantees can be derived when Assumption 1 is relaxed and replaced by suitable moment bounds on the rewards. We defer a detailed treatment of one such extension (in which a bound is assumed on the squared coefficient of variation of the maximum reward) to the appendix, Section A.5.

3. Polynomial-Time Algorithm and Its Analysis

In this section, we describe and analyze our main algorithm $\mathcal{A}$, which (approximately) implements the stopping policy described informally in Section 2.3, that is, stop the first time $Z_t^k\ge -\frac{U}{k-1}$, with $k=O\left(\frac{U}{ϵ}\right)$ and $Z_t^k$ approximated by Monte Carlo simulation (involving depth k nested conditional expectations). Because by (5) and the definition of MAR, $Z_t^k=Z_t-E[{\sum }_{j=1}^{k-1}{\max }_{i=1,\dots ,T}{Z}_i^j|{\mathcal{F}}_t]$, such a stopping time is equivalent to stopping the first time $Z_t-E[{\sum }_{j=1}^{k-1}{\max }_{i=1,\dots ,T}{Z}_i^j|{\mathcal{F}}_t]\ge -\frac{U}{k-1}$ (with $k=O\left(\frac{U}{ϵ}\right)$). Thus, $E[{\sum }_{j=1}^{k-1}{\max }_{i=1,\dots ,T}{Z}_i^j|{\mathcal{F}}_t]$ serves as an expansion-based approximation to MAR, with our stopping time the corresponding approximation to ${\tau }_{\infty }$.

We note that given such an algorithm, one can easily derive an algorithm for approximating the optimal value (by generating random sample paths, implementing $\mathcal{A}$ on them, and averaging/taking the median), which we later state as a corollary. We begin by formalizing a generative/simulation model for the OS problem.

Let us also formalize a computational model for analyzing algorithms. We adopt a model in line with Shapiro and Nemirovski (2005) and Swamy and Shmoys (2012), in which one is given access to an appropriate (simulation) oracle and can perform basic arithmetic operations on real numbers (regardless of the number of bits). Such a model falls within the framework of information complexity (Traub and Werschulz 1998), and we refer the reader to De Klerk (2008) for a more detailed comparison with the Turing-machine/bit complexity model.

Consistent with the assumptions in these works, we suppose that addition, subtraction, multiplication, division, maximum, and minimum of any two numbers can be done in one unit of time, as can the raising of one number to the power of another, regardless of the values of those numbers. We also suppose that the median of n numbers can be computed in n units of computational time, using, for example, the celebrated linear-time method for computing the median (Cormen et al. 2009), and ignoring the absolute constant governing the associated linear run time. We ignore all computational costs associated with reading, writing, and storing numbers in memory, as well as inputting numbers to functions. In general, our model allows for computation over real numbers. We also assume that drawing a sample from a standard Bernoulli, uniform ($\mathcal{U}[0,1]$), or normal ($\mathcal{N}(0,1)$) r.v.s. takes one unit of time. Most importantly, we assume calling $\mathcal{B}$ to generate a (conditioned) trajectory of $\mathbf{Y}$ takes C units of time, and evaluating the reward function $g_t(·)$ takes G units of time. We assume that both G and C scale at most polynomially with T and D.

3.1. Algorithm Description and Guarantees

Before stating algorithm $\mathcal{A}$, we first formalize a family of subroutines $\{{\mathcal{B}}^k,k\ge 1\}$ that can recursively (approximately) compute the value of $Z_t^k$ conditional on $Y_{[t]}=\gamma$ for any given input $\gamma \in {\mathbb{R}}^{D\times t}$ (which we denote by $Z_t^k(\gamma )$) to accuracy $ϵ$ with probability at least $1-\delta$ (with $ϵ,\delta$ inputs to ${\mathcal{B}}^k$). When $k=1$, ${\mathcal{B}}^1(t,\gamma )$ simply returns the exact reward value $g_t(\gamma )$. When $k\ge 2$, as $Z_t^k(\gamma )=Z_t^{k-1}(\gamma )-E[{\max }_{i=1,\dots ,T}{Z}_i^{k-1}|Y_{[t]}=\gamma ],$ the algorithm ${\mathcal{B}}^k$ approximates $Z_t^k(\gamma )$ by natural Monte Carlo simulation of the conditional expectation, with repeated recursive calls to ${\mathcal{B}}^{k-1}$. A detailed description of ${\mathcal{B}}^k$ can be found in Algorithm 1. Here, we recall that the “median trick” (as we formalize in Lemma A.2 of the appendix, Section A.4) is a direct corollary of Hoeffding’s inequality that is commonly used in speeding up estimation algorithms (by returning the median of independent estimations; see, e.g., Wang and Han 2015), and which goes back at least to the work of Nemirovski and Yudin (1983). By carefully adjusting (as a function of k) the parameters $ϵ$ and $\delta$ at different depths of recursion, we can ensure that ${\mathcal{B}}^k$ achieves the desired guarantees. We formalize the performance guarantee of ${\mathcal{B}}^k$ in the following result, deferring the proof to Section 3.2.

Theorem 3 demonstrates an interesting contrast between our expansion-based approach and standard methods such as DP and naive nested simulation. To achieve a similar performance guarantee, those methods typically require an exponential-in-T run time. In contrast, our approach yields polynomial-in-T algorithms that can achieve any fixed-precision guarantee, in the precise sense that once $ϵ$ is fixed, the exponent of T in the stated complexity bounds is a constant (depending on $ϵ$). However, let us point out that $\mathcal{A}$ is not a practically efficient algorithm, as is demonstrated in our numerical experiments (see the Online Companion). Our results thus serve as a proof of concept that a polynomial-time algorithm exists for complex OS problems. We leave the development of practical algorithms inspired by these theoretical insights as a direction for future research.

We end this section by noting that algorithm $\mathcal{A}$ also implies an efficient algorithm for approximating the optimal value with high probability (by generating random sample paths, implementing $\mathcal{A}$ on them, and averaging/taking the median). We provide a formal proof in the appendix, Section A.3.

3.2. Analysis

We analyze ${\mathcal{B}}^k$ and $\mathcal{A}$ in this section, thereby proving Lemma 6 and Theorem 3.