Multistage Conditional Compositional Optimization Buse Şen† , Yifan Hu‡ , and Daniel Kuhn†
arXiv:2604.14075v1 [math.OC] 15 Apr 2026
† EPFL, Switzerland
[email protected] [email protected] ‡ Rutgers University, USA [email protected]
April 16, 2026
Abstract We introduce Multistage Conditional Compositional Optimization (MCCO) as a new paradigm for decision-making under uncertainty that combines aspects of multistage stochastic programming and conditional stochastic optimization. MCCO minimizes a nest of conditional expectations and nonlinear cost functions. It has numerous applications and arises, for example, in optimal stopping, linear–quadratic regulator problems, distributionally robust contextual bandits, as well as in problems involving dynamic risk measures. The naı̈ve nested sampling approach for MCCO suffers from the curse of dimensionality familiar from scenario tree-based multistage stochastic programming, that is, its scenario complexity grows exponentially with the number of nests. We develop new multilevel Monte Carlo techniques for MCCO whose scenario complexity grows only polynomially with the desired accuracy.
1
Introduction
Multistage stochastic optimization constitutes a central and still largely unresolved challenge in the field of optimization under uncertainty. Efficient exact solution methods are out of reach even if the probability distribution of the random problem parameters is given in closed form. Indeed, even linear two-stage stochastic programs with independent and uniformly distributed random parameters are #P-hard [20, 34]. Consequently, multistage stochastic programs are often studied in a black-box framework, where the decisionmaker has no analytical distributional information and aims to construct an ϵ-optimal solution relying only on oracles that supply samples of the random problem parameters. In this framework, multistage stochastic programs are often addressed with the Sample Average Approximation (SAA), which builds a scenario tree using conditional sampling and solves a finite-dimensional approximate stochastic program on this empirically generated tree. However, the scenario complexity of SAA for obtaining an ϵ-optimal solution to a T -stage linear stochastic program is bounded above by O(ϵ−2T ) [56, 57], and this bound cannot be significantly improved by any SAA-type method even if the random parameters are serially independent [52, Proposition 2]. Hence, the size of this scenario tree grows exponentially with the number of stages. We are not aware of any other black-box method that overcomes this curse of dimensionality. With this paper we hope to lay the groundwork for a new generation of black-box methods for multistage stochastic optimization whose scenario complexity depends on the desired accuracy ϵ only through (1/ϵ)p , where the exponent p is independent of T . To this end, we study Multistage Conditional Compositional Optimization (MCCO) problems, which are reminiscent of multistage stochastic programs. In fact, we 1
will show that MCCO encapsulates various stochastic dynamic control problems with serially dependent disturbances as special cases. In addition, we will argue that MCCO is highly expressive and even includes problem instances that cannot be reformulated as any risk-neutral multistage stochastic program. The main contribution of this paper is to propose new multilevel Monte Carlo methods for MCCO problems with T ≥ 2 stages whose scenario complexity scales as O(ϵ−2 ) under standard regularity conditions. The MCCO problems to be investigated in this paper can be represented as min F (x), x∈X
(MCCO)
where X ⊆ Rd denotes the feasible set, and where the objective function F : Rd → R is defined as a nest of T conditional expectations and nonlinear integrands of the form " !# h i . (MCCE) F (x) := E f1 ξ1 , E1 f2 ξ2 , . . . ET −2 fT −1 (ξT −1 , ET −1 [fT (ξT , x)]) . . . Here, ξ1 , . . . , ξT represents a stochastic process, E[·] denotes unconditional expectation, and Et [·] denotes expectation conditional on the σ-algebra Ft generated by the time-t observation history ξ[t] := (ξ1 , ..., ξt ). The integrands ft : Rmt × Rdt → Rdt−1 with d0 = 1 and dT = d can be interpreted as uncertainty-affected cost functions. A key assumption throughout the paper is that the random vectors ξt are dependent and that the integrands ft are nonlinear. Hence, evaluating the objective function at a given x ∈ X is challenging, a problem which we henceforth refer to as Multistage Conditional Compositional Estimation (MCCE). MCCO reduces to classical (single-stage) stochastic optimization when T = 1, to conditional stochastic optimization when T = 2 [18, 29, 37, 39–41], to multilevel compositional stochastic optimization when the random vectors ξt are mutually independent [4, 14, 17, 21, 44, 54, 69–71] and to MCCE when the feasible set X is a singleton [51, 61, 72]. As we allow T to be any integer and the random vectors ξt to be dependent, MCCO strictly generalizes all of these more basic decision-making paradigms. Both MCCE and MCCO arise in many applications relevant for operations research, economics or engineering. For example, the problem of computing credit valuation adjustments can be naturally formulated as an instance of MCCE [28]. MCCE also emerges in optimal stopping problems, which encompass applications such as option pricing in finance [72] or inventory replenishment in operations management [50]. In addition, we will show that linear quadratic regulator problems and general finite-state stochastic control problems can be reformulated as instances of MCCE. Remarkably, these reformulations are possible even when the disturbances are serially dependent, as is the case when the underlying system is in a transient regime or is subject to distribution shifts. Finally, MCCO naturally arises in distributionally robust contextual bandits [59] or nested risk minimization [67]. All of these examples, which will be presented in detail in Section 2, underscore the strong ties between MCCE, MCCO and multistage stochastic optimization. However, MCCO is not subsumed by multistage stochastic optimization. Indeed, the infimum of any riskneutral multistage stochastic program is concave in the joint distribution of the random problem parameters, whereas F (x) clearly fails to be concave in the distribution of ξ[T ] for some choices of the integrands ft . Much like in multistage stochastic optimization, efficient exact solution methods for MCCE and MCCO are out of reach even if the probability distribution of the underlying stochastic process is given in closed form. Indeed, solving MCCE exactly is #P-hard even if T = 1, ξ1 follows the uniform distribution on the unit hypercube in Rm1 and f1 is a piecewise affine function with only two pieces [34, Corollary 1]. We thus adopt a black-box model, that is, we assume access to T oracles (“black boxes”) that one can use to generate any number of independent samples from the marginal distribution of ξ1 and the conditional distribution of ξt+1 given ξ[t] , t ∈ [T − 1], respectively. However, we do not presume to know these distributions. In this situation, we can resort to sampling methods for solving MCCE and MCCO approximately.
2
In this paper we will design and analyze function value estimators Fb(x) for F (x) that attain a specified mean squared error with a minimum number of scenarios (i.e., sample paths), and we will refer to this number as the estimator’s mean squared error-based scenario complexity. Given a function value estimator Fb(x), one can attempt to solve the following sample-based approximation of the intractable MCCO problem. min Fb(x) x∈X
(1)
We will show that the minimizers of problem (1) are near-optimal in (MCCO) with high probability. However, computing these minimizers can be challenging in the absence of restrictive convexity and monotonicity conditions. Thus, we will also design and analyze sample-efficient estimators for ∇F (x) that can be employed within stochastic gradient descent-type methods for solving problem (MCCO) to stationarity. The SAA estimator is arguably the simplest estimator for F (x). It is obtained by replacing all conditional expectations in (MCCE) with conditional sample averages. This amounts to approximating the stochastic process of the disturbances ξt , t ∈ [T ], by a scenario tree. The number of samples used to approximate Et [·] coincides with the scenario tree’s branching factor at stage t. We will show that the scenario complexity of an SAA estimator with optimally tuned branching factors amounts to O(ϵ−2T ) if all integrands ft are Lipschitz continuous and to O(ϵ−(T +1) ) if the integrands are additionally smooth. In both cases, the SAA estimator suffers from a curse of dimensionality. To design more efficient estimators, we leverage Multilevel Monte Carlo (MLMC) methods, which have originally been popularized in computational finance and uncertainty quantification [23, 24]. A basic MLMC estimator can be constructed as follows. For any level ℓ ∈ N0 , let Ybℓ denote the SAA estimator of F (x) with a uniform branching factor 2ℓ across all nodes of the scenario tree, and set Yb−1 = 0. By the law of large numbers, Ybℓ converges surely to F (x) as ℓ grows. Hence, F (x) P∞ almost b b admits the telescoping sum representation F (x) = E[ ℓ=0 (Yℓ − Yℓ−1 )]. If λ is a discrete random variable independent of all samples and if P(λ = ℓ) = q(ℓ) > 0 for all N0 , then the law of total expectation implies "∞ # # " X Ybℓ − Ybℓ−1 Ybλ − Ybλ−1 q(ℓ) F (x) = E =E . (2) q(ℓ) q(λ) ℓ=0
The inverse propensity weighted telescoping term Fb(x) = (Ybλ − Ybλ−1 )/q(λ) associated with level λ constitutes a randomized MLMC estimator in the spirit of [11, 53]. By construction, it satisfies E[Fb(x)] = F (x) and is thus unbiased. The choice of the distribution for λ is critical for controlling its computational cost. To see this, note that the SAA estimator Ybℓ requires 2ℓ(T −1) scenarios—one for each leaf node of the underlying scenario tree—implying that its computational cost grows quickly with ℓ. Drawing λ from a geometric distribution, for instance, ensures that small (inexpensive) levels ℓ are sampled with high probability, while large (expensive) levels are sampled rarely. This trick allows us to reduce the expected computational cost of Fb(x) while preserving its unbiasedness. To further ensure that the variance of Fb(x) remains bounded, we reuse all samples of the level ℓ − 1 scenario tree when constructing the level ℓ tree. Hence, Ybℓ and Ybℓ−1 are strongly positively correlated, and the variance of Fb(x) remains bounded thanks to a control variate effect. In this paper, we introduce new MLMC estimators for F (x) and ∇F (x) whose construction departs from that of the naı̈ve MLMC estimator described above in three key respects. First, we truncate the telescoping sum in (2) at a maximum level M in order to control the estimator’s variance [10]. Truncation introduces a bias, but the bias decays exponentially with M . Hence, selecting M becomes a matter of balancing bias and variance. Second, we adopt a stagewise recursive construction employing antithetic sampling [27], thereby improving the estimator’s computational cost. Third, we average multiple independent and identically distributed (i.i.d.) copies of the same estimator to reduce the variance. We will prove that the scenario complexity of the resulting MLMC estimator for F (x) amounts to O(log(ϵ−1 )2T −2 ϵ−2 ) if all integrands ft are Lipschitz continuous and to O(ϵ−2 ) if the integrands are additionally smooth. Under an additional Hölder continuity assumption, analogous results are obtained for our MLMC estimator of ∇F (x). The main contributions of this paper can be summarized as follows. 3
1. SAA value estimator: We provide a refined scenario complexity analysis of the SAA estimator for F (x). First, we prove that the mean squared error-based scenario complexity of this estimator scales as O(ϵ−2T ) and O(ϵ−(T +1) ) for nonsmooth and smooth integrands, respectively. Under a standard sub-Gaussianity assumption, we further show that, for nonsmooth and smooth integrands respectively, O(log(ϵ−1 )ϵ−2T ) and O(log(ϵ−1 )ϵ−(T +1) ) scenarios suffice to ensure that the estimator differs from F (x) at most by ϵ, uniformly across all x ∈ X and with high probability. All bounds are presented explicitly without any unknown constants. Thus, they provide finite-sample guarantees. 2. MLMC value estimator: We propose a new recursive truncated MLMC estimator for F (x) together with sharp1 scenario complexity guarantees. First, we prove that the mean squared error-based scenario complexity of this estimator scales as O(log(ϵ−1 )2T −2 ϵ−2 ) and O(ϵ−2 ) for nonsmooth and smooth integrands, respectively. Under a standard sub-Gaussianity assumption, we further show that, for nonsmooth and smooth integrands respectively, O(log(ϵ−1 )T ϵ−2 ) and O(log(ϵ−1 )ϵ−2 ) scenarios suffice to ensure that the estimator differs from F (x) at most by ϵ, uniformly across all x ∈ X and with high probability. Consequently, we provide the first scenario complexity results for MCCE with nonsmooth integrands that eliminate the usual exponential dependence on T in ϵ−1 . 3. MCCO paradigm: We introduce MCCO as a new paradigm for optimization under uncertainty and show that it unifies and generalizes numerous traditional decision-making models from the literature. 4. MLMC gradient estimator: We propose new recursive truncated as well as untruncated MLMC estimators for ∇F (x) together with sharp scenario complexity guarantees. Under the assumption that the integrands are smooth and have Hölder continuous Hessians, we prove that the mean squared errorbased scenario complexities of both estimators scale as O(ϵ−2 ). If the Hölder continuity assumption is relaxed, the scenario complexity increases to O(log(ϵ−1 )2T −2 ϵ−2 ) for the truncated and diverges for the untruncated estimator. We develop the first MLMC gradient estimators for MCCO with T ≥ 2. 5. SGD algorithm: We integrate our new MLMC gradient estimators into a stochastic gradient descent (SGD) algorithm for finding a stationary point of problem (MCCO). We prove that if all integrands are smooth and have Hölder continuous Hessians, then this algorithm needs O(ϵ−4 ) scenarios to reach an ϵ-stationary point both when equipped with the truncated and untruncated variant of the estimator. 6. Numerical experiments: We validate our theoretical results with numerical experiments revolving around Bermudan option pricing (for MCCE) and distributionally robust contextual bandits (for MCCO). These experiments confirm that the proposed MLMC estimators achieve the same mean square error as an SAA estimator using a significantly smaller number of scenarios. Literature Review MCCE with dependent disturbances has been widely studied in the literature. However, most works focus on problems with T = 2 stages, which arise in derivative pricing [13, 33] or risk management [12, 24, 26, 32]. For T = 2, it is well known that SAA estimators typically achieve a root mean squared error of ϵ using O(ϵ−3 ) scenarios, while MLMC estimators require only O(ϵ−2 ) scenarios [25]. The literature on MCCE with T > 2 stages is far more limited. It is shown in [51] that the scenario complexity of SAA estimators amounts to O(ϵ−2T ) for nonsmooth integrands and to O(ϵ−(T +1) ) for smooth ones. This work, however, provides only asymptotic rates without explicit constants. Moreover, it only studies mean squared errors but provides no uniform convergence analysis, which would be necessary to ensure that the minimizers of (1) approximate those of (MCCO) with high probability. To the best of our knowledge, only two papers develop MLMC estimators for MCCE with T > 2. The first focuses on 1
Some of these guarantees are sharp only up to poly-logarithmic factors in ϵ−1 .
4
optimal stopping problems and imposes restrictive assumptions on the integrands [72]. The second primarily considers smooth integrands. While it also discusses nonsmooth integrands, the corresponding MLMC estimators have infinite variance, making problem (1) numerically unstable [61]. No prior work develops MLMC estimators for MCCE with T > 2 stages and general nonsmooth integrands. To our best knowledge, all existing work on MCCO assumes that T = 2. In this case, MCCO reduces to conditional stochastic optimization, which arises in many applications. In machine learning, for instance, it appears in reinforcement learning [19, 39], model-agnostic meta-learning, robust supervised learning or invariant learning [39, 41]. In statistics and econometrics, applications include instrumental variable regression [8, 15, 29, 36, 49, 60] or Bayesian experimental design [5, 30, 31]. The simplest approach to MCCO is to approximate F (x) by an SAA estimator Fb(x) and to solve (1). For T = 2, the number of scenarios needed to guarantee that the resulting minimizer of (1) constitutes an ϵ-optimal solution to (MCCO) with high probability amounts to O(ϵ−4 ) if the integrands are nonsmooth and can be reduced to O(ϵ−3 ) if the integrands are smooth [39]. However, solving (1) may be computationally infeasible without restrictive convexity or monotonicity assumptions. To address this, subsequent work develops SAA estimators for ∇F (x) that can be employed within a stochastic gradient descent algorithm for solving problem (MCCO) to stationarity. For smooth integrands, this method requires O(ϵ−6 ) scenarios to find an ϵ-stationary point, or O(ϵ−5 ) with standard variance reduction techniques [40]. MLMC estimators for ∇F (x) introduced in [42] further reduce this to O(ϵ−4 ) scenarios. MCCO with T > 2 stages has not yet been investigated. The remainder of the paper is structured as follows. Section 2 illustrates motivating applications of MCCE and MCCO. Sections 3 and 4 then introduce and analyze SAA and MLMC estimators for F (x) and ∇F (x), respectively. Numerical results are presented in Section 5, and all proofs as well as some auxiliary technical results are relegated to the appendix. P Pn p 1/p . Notation The entrywise p-norm of a matrix A ∈ Rn×m is denoted by ∥A∥p := ( m j=1 i=1 |Aij | ) If p = 2, ∥A∥p reduces to the Frobenius norm. We use N0 and N to denote the sets of natural numbers with and without 0, respectively, and for any n ∈ N, we set [n] = {1, 2, . . . , n}. The transpose of the Jacobian of a differentiable function f : Rn → Rm is denoted by ∇f : Rn → Rn×m . If m = 1, ∇f thus reduces to the gradient of f . All (in)equalities involving random variables are understood to hold almost surely. The variance of a random variable ξ is denoted by V(ξ). We designate estimators by the superscript ‘b’ and define C(Fb) as the number of scenarios (sample paths) needed to construct an estimator Fb. The truncated geometric distribution with rate r ∈ (0, 1) and truncation point M ∈ N is denoted by Geo(r|M ). Thus, a random variable λ governed by Geo(r|M ) satisfies P(λ = ℓ) = r(1 − r)ℓ /(1 − (1 − r)M +1 ) for every ℓ = 0, . . . , M . If a(ϵ) and b(ϵ) are real-valued functions of a tolerance ϵ > 0, then a(ϵ) = O(b(ϵ)) means that there exist C > 0 and ϵ̄ > 0 such that a(ϵ) ≤ Cb(ϵ) for all ϵ ∈ (0, ϵ̄). Note that a(ϵ) and b(ϵ) could depend on other variables. However, O(·) captures only their asymptotic dependence on ϵ.
2
Motivating Applications
Instances of (MCCE) and (MCCO) with T > 2 stages arise in numerous applications. Examples include optimal stopping, linear-quadratic regulator and finite-state stochastic control problems for (MCCE) as well as policy learning in contextual bandits and risk measurement with nested risk measures for (MCCO). Optimal Stopping Suppose that an agent observes a sequence of random variables ξt , t ∈ [T ], which are revealed one at a time, and let gt : Rmt → R, t ∈ [T ], be Borel measurable payoff functions. At each time t ∈ [T ], the agent must decide whether to stop observing, in which case an immediate payoff gt (ξt ) is received, or to continue observing, in which case no immediate payoff is received. The process must stop at time T at the latest. If the goal is to maximize the expected payoff, then the optimal stopping time can be 5
found via dynamic programming [72]. To this end, use Ut (ξ[t] ) to denote the agent’s expected future payoff conditional on all observations up to time t. One readily verifies that UT (ξ[T ] ) = gT (ξT ) and Ut (ξ[t] ) = max gt (ξt ), Et Ut+1 (ξ[t+1] )
∀t ∈ [T − 1].
Unraveling this recursion yields h n oi . E [U1 (ξ1 )] = E max g1 (ξ1 ), E1 max g2 (ξ2 ), . . . , ET −1 [gT (ξT )] . . . Thus, the problem of computing the expected value of the optimal stopping time can be viewed as an instance of (MCCE) with fT (ξT , xT ) = gT (ξT ) and ft (ξt , xt ) = max{gt (ξt ), xt } for all t ∈ [T − 1]. Optimal stopping problems emerge in various applications across many domains [16]. For example, they arise in the pricing of American-style options [72] or in sequential hypothesis testing [47]. In addition, they also arise in supply chain management [63] to determine the timing of inventory replenishment [50], and in game theory to model the reasoning about the probabilistic strategies of players [51]. Consider the finite-horizon stochastic linear-quadratic regulator problem # "T −1 X ⊤ ⊤ s⊤ min E t Qst + at Rat + sT PT sT a,s t=1 (3) s.t. st+1 = Ast + Bat + ξt+1 ∀t ∈ [T − 1] at is Ft -measurable ∀t ∈ [T − 1],
Linear-Quadratic Regulator
where st ∈ Rm and at ∈ Rn denote the state and the control action at time t. We assume that the system matrices A ∈ Rm×m and B ∈ Rm×n as well as the symmetric positive semidefinite cost matrices Q ∈ Rm×m , R ∈ Rn×n and PT ∈ Rm×m are deterministic and that R is invertible. We also assume that s1 = ξ1 and that {ξt }t∈N is a (possibly serially dependent) stochastic process with finite second moments taking values in Rm . Control problems of the form (3) are ubiquitous in economics [35], engineering [9] or neuroscience [62], among many others. Traditionally, the stochastic disturbances are assumed to be i.i.d. However, in an era where control systems must navigate continuously shifting environments and interact with adaptive, datadriven agents, models that abandon the fragile i.i.d. noise assumption are essential for ensuring robustness under real-world nonstationarity and persistent distribution shift. The Q-function at time t associated with problem (3) is defined via the Bellman-type recursion ⊤ ⊤ ′ Qt (s, a) = s Qs + a Ra + Et min Qt+1 (As + Ba + ξt+1 , a ) (4) ′ n a ∈R
with terminal value QT (s, a) = s⊤ PT s (which is deterministic and independent of a); see also [9, § 4.1]. Due to the serial dependencies of the noise process, Qt (s, a) is an Ft -measurable random object for any fixed (s, a). But we notationally suppress the dependence of Qt (s, a) on ξ[t] to avoid clutter. The optimal value of (3) can be expressed in terms of the time-1 Q-function as E minn Q1 (ξ1 , a) . (5) a∈R
The following lemma shows that the Q-functions are quadratic in the states and actions with random coefficients. To keep the paper self-contained, we sketch the proof of this elementary result in Appendix B.
6
Lemma 2.1. For each t ∈ [T ], the Q-function is quadratic in (s, a) and can be represented as ⊤ s ⊤ Qss Qsa s s b s t t Qt (s, a) = + 2 at + ct ⊤ Qaa a a b a (Qsa ) t t t
P-a.s.,
(6a)
where the Ft -measurable coefficient matrices and vectors are defined through the backward recursions ⊤ Qss t = Et Q + A Pt+1 A ⊤ Qsa t = Et A Pt+1 B ⊤ Qaa t = Et R + B Pt+1 B (6b) bst = Et A⊤ Pt+1 ξt+1 + A⊤ gt+1 bat = Et B ⊤ Pt+1 ξt+1 + B ⊤ gt+1 ⊤ ⊤ ct = Et ξt+1 Pt+1 ξt+1 + 2 gt+1 ξt+1 + dt+1 for all t ∈ [T − 1] initialized by Qss T = PT ,
Qsa T = 0,
Qaa T = 0,
bsT = 0,
cT = 0.
(6c)
−1 a dt = ct − (bat )⊤ (Qaa t ) bt
(6d)
baT = 0
and
The underlying auxiliary variables are defined by sa aa −1 sa ⊤ Pt = Qss t − Qt (Qt ) (Qt ) ,
aa −1 a gt = bst − Qsa t (Qt ) bt
and
for all t ∈ [T − 1], while PT is the terminal state cost matrix, gT = 0 and dT = 0. For all t ∈ [T ] we define an aggregate coefficient vector encoding the time-t Q-function as sa aa s a d xt = Qss t , Qt , Qt , bt , bt , ct ∈ R , where the components are the Ft -measurable matrices and vectors appearing in (6a), and where d = m2 + mn + n2 + m + n + 1. Also, for all t = 2, . . . , T we define ft : Rm × Rd → Rd through ft (ξt , xt ) = Q+A⊤ Pt A, A⊤ Pt B, R+B ⊤ Pt B, A⊤ Pt ξt +A⊤ gt , B ⊤ Pt ξt +B ⊤ gt , ξt⊤ Pt ξt +2gt⊤ ξt +dt , where Pt , gt and dt are defined as in (6d) and thus constitute rational functions of xt . In addition, we define f1 : Rm × Rd → R through ⊤ sa aa −1 sa ⊤ aa −1 a −1 a f1 (ξ1 , x1 ) = ξ1⊤ Qss ξ1 + 2 bs1 − Qsa ξ1 + c1 − (ba1 )⊤ (Qaa 1 − Q1 (Q1 ) (Q1 ) 1 (Q1 ) b1 1 ) b1 . The backward recursions (6b) for the Q-function coefficients now take the compact form xt = Et ft+1 (ξt+1 , xt+1 ) ∀t ∈ [T − 1] with terminal condition xT = (PT , 0, 0, 0, 0, 0). Using similar arguments as in the proof of Lemma 2.1, one readily verifies that the optimal value (5) of the control problem (3) can be represented compactly as E[f1 (ξ1 , x1 )]. In view of the above recursion, this implies that the problem of computing the optimal value of (3) is equivalent to solving problem (MCCE), where x is identified with the deterministic vector xT . A similar reasoning as shown here can be used to reduce linear exponential quadratic regulator and risk-sensitive control problems [68] to instances of (MCCE). We omit the derivations for brevity.
7
Stochastic Control with Finite State and Action Spaces We use the same symbols as in the last example to denote states, actions and disturbances and consider the finite-horizon stochastic control problem "T −1 # X min E ct (st , at ) + cT (sT ) a,s
s.t.
t=1
at ∈ At (ξt ) ∀t ∈ [T − 1] st+1 = ψ(st , at , ξt+1 ) ∀t ∈ [T − 1] at is Ft -measurable ∀t ∈ [T − 1].
(7)
All states and actions take values in finite sets S and A, respectively, while {ξt }t∈N represents a generic (possibly serially dependent) stochastic process of disturbances taking values in Rm . The system evolves according to a generic state transition function ψ : S × A × Rm → S, and at each time t, the controller must select an action from a generic random and time-varying feasible set At (ξt ) ⊆ A. The stage-costs ct : S × A → R, t ∈ [T − 1], and the terminal cost cT : S → R are also captured by generic functions. If the disturbances are serially independent, then problem (7) reduces to a Markov decision process, a problem class that occupies central stage in operations research and machine learning [9]. Allowing the disturbances to display serial correlations enlarges this problem class significantly. The time-t Q-function of problem (7) is defined by the Bellman-type recursion Qt (s, a) = ct (s, a) + Et min Qt+1 ψ(s, a′ , ξt+1 ), a′ (8) a′ ∈At+1 (ξt+1 )
with terminal condition QT (s, a) = cT (s). As in the previous example, Qt (s, a) is an Ft -measurable random object for any fixed (s, a), but we suppress its dependence on ξ[t] to avoid clutter. The optimal value of problem (7) can be expressed in terms of the time-1 Q-function as E min Q1 (ξ1 , a) . (9) a∈A1 (ξ1 )
For all t ∈ [T ] we define xt ∈ Rd with d = |S| × |A| as the vectorized time-t Q-function with components (xt )(s,a) = Qt (s, a), (s, a) ∈ S × A. Also, for all t = 2, . . . , T , we define ft : Rm × Rd → Rd through ft (ξt , xt ) (s,a) = ct−1 (s, a) + min xt ψ(s, a, ξt ), a′ ∀(s, a) ∈ S × A. a′ ∈At (ξt )
In addition, we define f1 : Rm × Rd → R through f1 (ξ1 , x1 ) = mina′ ∈A1 (ξ1 ) x1 (s, a′ ), where we use again the convention that s = ξ1 . With this notation, the Bellman recursions (8) can be written compactly as xt = Et ft+1 (ξt+1 , xt+1 ) ∀t ∈ [T − 1], with terminal condition (xT )(s,a) = cT (s) for every s ∈ S and a ∈ A. The optimal value (9) of problem (7) thus coincides with E[f1 (ξ1 , x1 )], and computing this value is equivalent to solving problem (MCCE). Distributionally Robust Policy Learning for Offline Contextual Bandits The contextual bandit problem is a sequential decision problem under uncertainty. In each round, an agent first observes a context c ∈ C, then selects an action a ∈ A (metaphorically thought of as an arm of a slot machine) and finally incurs a random cost ya ∈ Y associated with the chosen action. The distribution Qya |c of ya conditional on c is unknown to the agent. However, the pairs (c, y) with y = (ya )a∈A are known to be independent and identically distributed across rounds. The crux of multi-armed bandit problems is that, in each round, the agent observes only the cost of the chosen action but not those of the other actions. Thus, each action a ∈ A must be selected repeatedly for each given c in order to learn the distribution Qya |c . 8
Consider now a parametric family of policies πθ , θ ∈ Θ, which map any context c ∈ C to a probability distribution over the actions a ∈ A. We will henceforth assume that A is finite. Thus, given a context c, the policy πθ selects action a with probability πθ (a|c). The agent aims to find θ ∈ Θ such that πθ minimizes the average cost. If the joint distribution Q(c,y) of c and y was known, this could be achieved by solving Z Z X Z X ya dQy|c (y) dQc (c). min ya πθ (a|c) dQ(c,y) (c, y) = min πθ (a|c) θ∈Θ
θ∈Θ
C×Y a∈A
C a∈A
Y
In off-policy learning, however, Q(c,y) is unknown, and the agent has only access to context, action and cost data collected under a fixed behavioral policy. Any estimator P(c,y) constructed from this data invariably differs from Q(c,y) . To combat the detrimental effects of estimation errors, it has been proposed to solve the following distributionally robust variant of the off-policy learning problem [59]. Z X Z min sup ya dQy|c (y) dQc (c) (10) πθ (a|c) sup θ∈Θ Qc ∈P(C): d(Qc ,Pc )≤rc
C a∈A
Qy|c ∈P(Y): d(Qy|c ,Py|c )≤ry
Y
Here, we use the 2-Wasserstein metric d(·, ·) to measure distances between probability distributions. Note that the distributionally robust off-policy learning problem replaces the expectations with respect to the unknown distributions Qc and Qy|c by worst-case expectations over all distributions in an rc -neighborhood around Pc and an ry -neighborhood around Py|c , respectively. By [46, Proposition 6.20], the maximization R problem over Qy|c in (10) can be solved in closed form, and its supremum equals Y ya dPy|c (y) + ry . In addition, by [46, Theorem 4.18], the maximization problem over Qc can be dualized and thus be converted to a minimization problem over a scalar dual decision variable λ. Hence, problem (10) can be recast as ! Z Z X min sup πθ (a|u) ya dPy|c=u (y) + ry + rc2 λ − λ∥u − c′ ∥22 dPc (c′ ). (11) θ∈Θ,λ≥0
C u∈C
Y
a∈A
Note that (11) involves a maximization problem over all contexts u ∈ C, which may be difficult to solve because the expected reward of action a conditional on c = u is a measurable function of u and because maximizing a generic measurable function is hard. As a remedy, the supremum over all u ∈ C can be approximated by a softmax function involving the uniform distribution Pu on C. Thus, (11) simplifies to 1 log C µ
Z min
θ∈Θ,λ≥0
Z exp µ C
X a∈A
Z πθ (a|u)
ya dPy|c=u (y) Y
!! + ry + rc2 λ − λ∥u − c′ ∥22
!
(12)
dPu (u) dPc (c′ ),
where the parameter µ > 0 controls the smoothness of the softmax function. For further details see [59]. Problem (12) is readily recognized as an instance of (MCCO) with T = 3, x = (θ, λ), X = Θ×R+ , ξ1 = c′ , ξ2 = u and ξ3 = (c′ , u, y). The integrands are given by f1 (ξ1 , x1 ) = µ1 log(x1 ), f2 (ξ2 , x2 ) = exp(µx2 ) and f3 (ξ3 , x) =
X a∈A
πθ (a|u) ya + ry + rc2 λ − λ∥u − c′ ∥22 .
Variants of the distributionally robust off-policy learning problem (10) use the 2-Wasserstein distance with a Kullback-Leibler regularization term [3] or the Sinkhorn divergence [65] instead of the plain 2-Wasserstein distance to measure the discrepancy of different distributions. These alternative formulations also give rise to instances of MCCO with T = 3. For further details we refer to [65, Appendix D]. 9
Nested Risk Measures Consider a stream of random costs ξt , for t ∈ [T ], that are p-integrable for some exponent p ≥ 1 (i.e., E[∥ξt ∥p2 ] < ∞ for all t ∈ [T ]). In this situation, it makes sense to assign the total cost PT t=1 ξt a risk that depends on the temporal structure of the underlying stochastic process. This can be done by using a nested dynamic risk measure ρ = ρ0 ◦ · · · ◦ ρT −1 , where ρt : Lpt+1 → Lpt is a conditional risk mapping and Lpt represents the family of all p-integrable functions of the history ξ[t] of costs up to time t. The defining properties of a conditional risk mapping are translation invariance, convexity and monotonicity; see [55] for details. By exploiting these properties, the risk of the total cost can then be represented as i hP T ρ t=1 ξt = ρ0 [ξ1 + ρ1 [ξ2 + · · · + ρT −1 [ξT ] · · · ]] . A widely studied choice is to set ρt to a conditional entropic risk measure, defined as ρt (ξt+1 ) =
1 log (Et [exp(−µt ξt+1 )]) , µt
where µt > 0 represents the risk aversion level at time t. Notably, the risk aversion level may vary over time. Such scenarios arise, for instance, in pension planning, where the goal is to assess the risk of a pension contract. In this context, tolerance for uncertainty in cashP flows typically decreases over time as retirement approaches. It is now easy to see that evaluating exp(µ0 ρ( Tt=1 ξt )) reduces to an instance of (MCCE), with ft (ξt , xt ) = exp(−µt−1 (ξt + µ−1 t log(xt ))) for t ∈ [T − 1] and fT (ξT , x) = exp(−µT −1 ξT ). Similarly, optimization problems involving nested risk measures can often be reformulated as instances of (MCCO). Such problems arise, for example, in Bayesian distributionally robust optimization [58] and in the analysis of Bayesian risk Markov decision processes [48], which emerge in betting or control problems [48, 67]. Solving Bayesian risk Markov decision processes presents significant challenges due to the complexity of computing nested risk functionals [48]. To address these difficulties, approximation techniques are commonly used. However, traditional estimation and optimization methods can be computationally expensive. This paper proposes cost-efficient methods to mitigate this issue and resolves an open question raised in [48].
3
Numerical Solution of MCCE Problems
We begin by presenting the fundamental assumptions that underpin this paper. Assumption 3.1 (Basic Regularity Conditions). (i) The integrand ft (ξt , xt ) is Borel-measurable in ξt , and E[∥ft (ξt , 0)∥2 ] < ∞ for every t ∈ [T ]. (ii) The feasible set X is convex and compact. Thus, it has a finite diameter DX . Neither F nor the underlying integrands ft , t ∈ [T ], are assumed to display any convexity properties. However, we assume these functions to be Lipschitz continuous, which is common practice in sample complexity analysis (see e.g., [57]). Sometimes, we additionally assume the integrands to be smooth. Assumption 3.2 (Lipschitz Continuity). ft (ξt , xt ) is Lt -Lipschitz continuous in xt uniformly for all ξt , i.e., ∥ft (ξt , xt ) − ft (ξt , yt )∥2 ≤ Lt ∥xt − yt ∥2
∀xt , yt ∈ Rdt , ∀ξt ∈ Rmt .
Assumption 3.3 (Smoothness). ft (ξt , xt ) is St -smooth in xt for across all ξt , i.e., ∥∇xt ft (ξt , xt ) − ∇xt ft (ξt , yt )∥2 ≤ St ∥xt − yt ∥2
10
∀xt , yt ∈ Rdt , ∀ξt ∈ Rmt .
Note that if Assumption 3.3 holds, then the fundamental theorem of calculus implies that ft (ξt , xt ) − ft (ξt , yt ) − ∇xt ft (ξt , yt )⊤ (xt − yt )
2
≤
St ∥xt − yt ∥22 2
∀xt , yt ∈ Rdt , ∀ξt ∈ Rmt . (13)
Assumptions 3.1 and 3.2 ensure that (MCCE) and (MCCO) are well-defined. Indeed, they imply that all conditional expectations in the definition of F (x) exist and are finite, which is an immediate consequence of Lemma C.2 in the appendix. Thus, problem (MCCE) is well-defined. Q In addition, Assumptions 3.1 and 3.2 imply that F (x) is Lipschitz continuous with Lipschitz constant Tt=1 Lt . This is an immediate corollary of Lemma C.3 in the appendix. By Assumption 3.1 (ii), we can thus invoke Weierstrass’ extreme value theorem to conclude that the minimum in (MCCO) is attained. Hence, problem (MCCO) is well-defined. Problem (MCCE) is challenging due to the inherent difficulty of high-dimensional integration. Consequently, attempting to compute F (x) exactly is overly ambitious. In this section we thus pursue the more modest but achievable goal of approximating F (x) with an estimator Fb(x) constructed from samples. Specifically, we will construct an SAA estimator in Section 3.1 and an MLMC estimator in Section 3.2. We will show that both estimators converge to F (x) in mean squared error and in probability, uniformly over all x ∈ X . Moreover, we will demonstrate that the MLMC estimator outperforms the SAA estimator in terms of both mean squared error-based and high-probability scenario complexity.
3.1
SAA Estimator
Perhaps the most straightforward estimator for F (x) is the SAA estimator defined below. Definition 3.4 (SAA Estimator). We construct a forest of n1 scenario trees as follows. The root nodes of the trees are indexed by i1 ∈ [n1 ]. The children of i1 are indexed by i[2] = (i1 , i2 ) for i2 ∈ [n2 ], and the children of i[2] are indexed by i[3] = (i[2] , i3 ) for i3 ∈ [n3 ], etc. Thus, the multi-index i[t] = (i1 , i2 , . . . , it ) of a stage-t node encodes the tree that accommodates the node i[t] as well as the path from the root node in that tree to the node i[t] ; see Figure 1. The SAA estimator for F (x) is then given by 1 Fb(x) := n1
n1 X
(14a)
Fbi1 (x),
i1 =1
where the term corresponding to the i1 -th scenario tree is set to nT n2 i X X 1 1 i fT ξT[T ] , x . . . f2 ξ2[2] , . . . , Fbi1 (x) := f1 ξ1i1 , n2 nT iT =1
i2 =1
∀i1 ∈ [n1 ],
(14b)
i
and where the sample ξt[t] associated with node i[t] is constructed by the forward recursion i[t−1] i i.i.d. ξt[t] ∼ Pξt · ξ1 = ξ1i1 , . . . , ξt−1 = ξt−1
∀i1 ∈ [n1 ], . . . , ∀it ∈ [nt ], ∀t ∈ [T ].
We introduce some key notations used throughout the paper. Given any t ∈ N0 with t + 1 ∈ [T ] and any b t [ψ(ξ[t+1] )] as shorthand for the sample average Borel-measurable function ψ(ξ[t+1] ), we henceforth use E i of ψ(ξ[t+1] ) over nt+1 samples ξt+1 , i ∈ [nt+1 ], drawn independently from the marginal distribution of ξ1 (if t = 0) or the conditional distribution of ξt+1 given ξ[t] (if t > 0). That is, we set b t [ψ(ξ[t+1] )] := E
1 nt+1 11
nt+1
X i=1
i ψ(ξ[t] , ξt+1 ).
(15)
b as shorthand for E b 0 [·]. Using these conventions, the SAA estimator of DefiBelow, we sometimes use E[·] nition 3.4 can be recast more compactly as h i b f1 ξ1 , E b 1 f2 ξ2 , . . . E b T −1 [fT (ξT , x)] . . . Fb(x) = E .
ξ1i1 ...
(i ,1)
ξ2 1 (i ,1,1) . . .
ξ3 1
(i ,1,i3 ) . . .
ξ3 1
(i ,1,n3 )
ξ3 1
...
(i ,i )
ξ2 1 2 (i ,i ,1) . . .
ξ3 1 2
ξ3 1 2 3 . . . ξ3 1 2 (i ,i ,i )
(i ,i ,n3 )
(i ,n2 )
ξ2 1
(i ,n2 ,1) . . . (i1 ,n2 ,i3 ). . . (i1 ,n2 ,n3 ) ξ3 ξ3
ξ3 1
Figure 1: Visualization of the i1 -th scenario tree underlying the SAA estimator when T = 3. Q Note that the SAA estimator of Definition 3.4 requires C(Fb(x)) = Tt=1 nt scenarios. In the following we will prove that, as the sample sizes nt , t ∈ [T ], tend to infinity, the estimator Fb(x) converges in mean squared error and in probability to F (x) uniformly across all x ∈ X . To this end, we first rewrite the mean squared error of Fb(x) with respect to F (x) as the sum of the squared bias and the variance of Fb(x). h 2 2 i MSE(Fb(x)) = E Fb(x) − F (x) = E Fb(x) − F (x) + V Fb(x)
(16)
This formula shows that, for any fixed the mean squared error of Fb(x) below ϵ2 by √ ϵ > 0, we can drive ensuring that |E[Fb(x)] − F (x)| ≤ ϵ/ 2 and V(Fb(x)) ≤ ϵ2 /2. In the following, we construct bounds on the bias and the variance of the SAA estimator and analyze their dependence on the sample sizes nt , t ∈ [T ]. These bounds are obtained under the assumption that all integrands have finite variance. Assumption 3.5. Each integrand ft (ξt , xt ) has a finite conditional variance in the sense that i h 2 σt2 := sup ess sup Et−1 ft (ξt , xt ) − Et−1 [ft (ξt , xt )] 2 < ∞ ∀t ∈ [T ]. xt ∈Rdt
We are now ready to derive an upper bound on the bias of the Q SAA estimator. This upper bound is expressed in terms of product Lipschitz constants of the form L[t] := ts=1 Ls for t ∈ [T ] and L[0] = 1. Lemma 3.6. If Assumptions 3.1, 3.2 and 3.5 hold, then the bias of the SAA estimator satisfies P T −1 L√[t] σt+1 if Assumption 3.3 does not hold, t=1 n b E F (x) − F (x) ≤ PT −1 L t+1St σ2 [t−1] t+1 if Assumption 3.3 holds. t=1 2nt+1 Lemma 3.6 shows that the bias of the SAA estimator decays faster with the sample sizes if the integrands are smooth (i.e., if Assumption 3.3 holds). If T = 2, then Lemma 3.6 recovers the bounds of [39, Lemma 3.1] for conditional stochastic optimization. Next, we bound the variance of the SAA estimator. Lemma 3.7. If Assumptions 3.1, 3.2 and 3.5 hold, then the variance of the SAA estimator satisfies V(Fb(x)) ≤
T −1 2 1 X L[t] σt+1 + σ1 . √ n1 nt+1 t=1
12
Remark 3.8. The variance bound of Lemma 3.7 can be strengthened when the integrands are smooth. However, the improved bound is cumbersome and does not decay faster with the sample sizes. To keep our analysis transparent, we therefore use the simpler bound of Lemma 3.7 even if the integrands are smooth. Combining Lemmas 3.6 and 3.7 with (16) shows that the mean squared error of Fb (x) satisfies P 2 T −1 L[t] σt+1 2 [t] σt+1 1 PT −1 L√ √ + + σ if Assumption 3.3 does not hold, 1 t=1 t=1 n n n t+1 1 t+1 MSE(Fb(x)) ≤ 2 P P 2 T −1 L[t−1] St σt+1 2 T −1 L[t] σt+1 √ + n11 if Assumption 3.3 holds. t=1 t=1 2nt+1 nt+1 + σ1 We define the mean squared error-based scenario complexity of any estimator for F (x) as the total number of scenarios needed in order to ensure a root mean squared error of at most ϵ uniformly across all x ∈ X . The above bound enables us to characterize the mean squared error-based scenario Q complexity of the SAA estimator. To this end, recall that the construction of Fb(x) requires a total of Tt=1 nt scenarios. The next theorem shows that setting the sample sizes to nt = O(1/ϵ2 ) for nonsmooth integrands and nt = O(1/ϵ) √ for smooth integrands, for all t = 2, . . . , T , ensures that the bias is bounded by ϵ/ 2. Moreover, choosing √ n1 = O(1/ϵ2 ) guarantees that the standard deviation is bounded by ϵ/ 2. This implies that the mean squared error-based scenario complexity of the SAA estimator is at most O(ϵ−2T ) for nonsmooth integrands and O(ϵ−(T +1) ) for smooth integrands. The proof of this theorem is straightforward and therefore omitted. Theorem 3.9 (Mean Squared Error-Based Scenario Complexity √ of the SAA Estimator). Suppose that Assumptions 3.1, 3.2 and 3.5 hold. Assume also that n1 = ⌈1 + 2 2σ1 /ϵ + 2σ12 /ϵ2 ⌉ and that l √ m 2L[t−1] σt (T −1) 2 if Assumption 3.3 does not hold, ϵ nt = l √ m 2 2L[t−2] St−1 σt (T −1) if Assumption 3.3 holds, 2ϵ for all t = 2, . . . , T . Then, the SAA estimator achieves a root mean squared error of at most ϵ uniformly across all x ∈ X . Thus, the mean squared error-based scenario complexity of the SAA estimator amounts at most to O(ϵ−2T ) if Assumption 3.3 does not hold and to O(ϵ−(T +1) ) if Assumption 3.3 holds. Theorem 3.9 shows that although smoothness improves the exponent in the accuracy parameter ϵ, the scenario complexity of the SAA estimator still exhibits an exponential dependence on the horizon T in both the smooth and nonsmooth cases. This highlights a fundamental limitation of nested estimation. Achieving high accuracy becomes increasingly sample-intensive as T grows. It is important to note that these guarantees hold only in expectation; that is, the root mean squared error of the SAA estimator is at most ϵ on average over repeated independent simulation runs. To obtain stronger guarantees that hold for a single run with high confidence, we next turn to high-probability bounds. These results are based on concentration inequalities and therefore require an additional standard sub-Gaussianity assumption. Definition 3.10 (Sub-Gaussian Random Vector). An m-dimensional random vector ξ with finite mean is called sub-Gaussian with variance proxy ζ 2 > 0 if E exp(λ⊤ (ξ − E[ξ])) ≤ exp ∥λ∥22 ζ 2 /2 ∀λ ∈ Rm . Assumption 3.11 (Sub-Gaussianity of Fbi1 (x)). For each i1 ∈ [n1 ], the estimator Fbi1 (x) corresponding to the i1 -th scenario tree in Definition 3.4 is sub-Gaussian with variance proxy ζ 2 > 0 for any fixed x ∈ X . Assumption 3.11 is hard to verify in general, but Remark 3.12 below outlines conditions when it holds. Remark 3.12. Assumption 3.11 holds under either of the following easily verifiable conditions. 13
(i) If ξt is sub-Gaussian with variance proxy ζt2 and ft (ξt , xt ) is jointly Lt -Lipschitz continuous in ξt and xt for all t ∈ [T ], then Assumption 3.5 is automatically satisfied [64, Exercise 2.5]. In addition, one readily verifies that Fbi1 (x) is sub-Gaussian. Indeed, Fbi1 (x) can be constructed by iteratively averaging and concatenating sub-Gaussian random variables and transforming them under Lipschitz continuous functions. By Lemma C.5 in the appendix, all of these operations preserve sub-Gaussianity. (ii) If there exists a constant M > 0 such that |f1 (ξ1 , x1 )| ≤ M for all ξ1 ∈ Rm1 and x1 ∈ Rd1 , then [64, Example 2.4] readily implies that Fbi1 (x) is sub-Gaussian with variance proxy ζ 2 = 4M 2 . We now establish the uniform convergence of the SAA estimator to the true objective function. Lemma 3.13. Suppose that Assumptions 3.1, 3.2, 3.5 and 3.11 hold and that ϵ > 0 is a prescribed tolerance. Assume also that the stagewise sample sizes satisfy l 2 m 4L[t−1] σt (T −1) if Assumption 3.3 does not hold, ϵ nt = l m 2 2L[t−2] St−1 σt (T −1) if Assumption 3.3 holds, ϵ for all t = 2, . . . , T . Then, we have d n1 ϵ2 b . P sup |F (x) − F (x)| > ϵ ≤ 2 8L[T ] DX /ϵ + 1 exp − 32ζ 2 x∈X Recall that ζ 2 is the variance proxy of Fbi1 (x), d and DX are the dimension and the diameter of X , respectively, and L[T ] is the Lipschitz constant of F (x). Lemma 3.13 establishes a high probability bound on the worst possible difference between the SAA estimator and the true objective function F (x) uniformly across all x ∈ X . It asserts that the probability that this difference exceeds any given threshold ϵ > 0 decays exponentially fast with ϵ and with the sample size n1 . The uniform convergence result of Lemma 3.13 not only characterizes the number of scenarios needed to approximately solve the estimation problem (MCCE), which is the central focus of this section, but also inspires a straightforward method to approximately solve the optimization problem (MCCO). Indeed, instead of minimizing the function F (x), whose values and gradients are difficult to compute exactly, one may minimize the SAA estimator Fb(x), that is, one may solve an instance of (1). Below, we will use Lemma 3.13 to show that the exact minimizers of the SAA problem (1) yield near-optimal solutions to (MCCO) with high probability. The practicality of this approach ultimately depends on the tractability of (1). Remark 3.14. Global optimization of the SAA estimator Fb(x) is tractable in settings where Fb(x) is convex and admits an efficient separation oracle. A sufficient condition for convexity of both F (x) and Fb(x) is that the integrand ft (ξt , xt ) is convex in xt for every t ∈ [T ] and non-decreasing in xt for every t ∈ [T −1]. These conditions ensure that the global minimum of (1) can be found efficiently. However, they are restrictive. We now show that the minimizers of (1) are near-optimal in (MCCO) with high probability. Corollary 3.15. Assume that all conditions of Lemma 3.13 hold for some ϵ > 0, and let x∗ and x b∗ be minimizers of (MCCO) and the SAA problem (1), respectively. Then, we have d n1 ϵ2 ∗ ∗ P F (b x ) − F (x ) > ϵ ≤ 4 8L[T ] DX /ϵ + 1 exp − . 128ζ 2 We define the high-probability scenario complexity of a method for solving (MCCO) as the total number of scenarios that this method needs to identify an ϵ-optimal solution of (MCCO) with probability 1 − β. 14
Theorem 3.16 (High-Probability Scenario Complexity of the SAA Estimator). Suppose that all conditions of Lemma 3.13 hold for some ϵ > 0 and that 128ζ 2 d log 8L[T ] DX /ϵ + 1 + log (4/β) n1 = ϵ2 for some β ∈ (0, 1). Then any solution of the SAA problem (1) constitutes an ϵ-optimal solution of (MCCO) with probability at least 1 − β. Thus, the high-probability scenario complexity of the SAA estimator equals O(log(ϵ−1 )ϵ−2T ) if Assumption 3.3 does not hold and O(log(ϵ−1 )ϵ−(T +1) ) if Assumption 3.3 holds. Theorem 3.16 generalizes [39, Corollary 4.2] to general T > 2. It shows that, like the mean squared error-based scenario complexity, the high-probability scenario complexity of the SAA method grows exponentially with T . This aligns with the observation in multistage stochastic programming that SAA suffers from a curse of dimensionality [52, 56]. Note also that the exponential growth rate is reduced by a factor of two if the integrands are smooth in the sense of Assumption 3.3. In contrast to Theorem 3.9, the bound on n1 depends now on the dimension d of the decision variable and the variance proxy ζ 2 of the tree-wise SAA estimator. This captures the cost of enhancing guarantees from convergence in expectation to convergence with high probability. If the SAA estimator Fb(x) further satisfies the Hölderian error bound condition proposed in [39, Assumption 4.1], then the bound on n1 in Theorem 3.16 can be made independent of the dimension d. In the special case where T = 2, the scenario complexity bounds of Theorem 3.16 reduce to the known bounds for conditional stochastic optimization reported in [39]. These bounds scale quartically with the inverse accuracy 1/ϵ for nonsmooth integrands and cubically for smooth integrands, demonstrating the efficiency of the SAA method when T is small. However, for larger values of T , the method becomes increasingly inefficient, underscoring the need for more advanced and scalable approaches to solve (MCCO).
3.2
MLMC Estimator
We now develop a recursive MLMC estimator for the nested expectation problem (MCCE). Our approach generalizes the unbiased MLMC method of [11] to arbitrary depths T ≥ 2, but introduces a controlled bias to simultaneously guarantee finite variance and computational cost under weaker conditions (e.g., Lipschitz continuity). This contrasts with the unbiased recursive MLMC estimator of [61], which requires smoothness to ensure a finite variance. To this end, we need additional notation and definitions. For any t + 1 ∈ [T ] and Borel-measurable function ψ(ξ[t+1] ), we continue to use the empirical conditional expectation operab t [ψ(ξ[t+1] )] with nt+1 samples defined in (15). We further introduce empirical conditional expectation tor E operators that run over 2ℓ samples as well as over the even or odd indices in [2ℓ ]. In other words, we define 2ℓ
b ℓ [ψ(ξ[t+1] )] E t
:=
b ℓ,o [ψ(ξ[t+1] )] := E t b ℓ,e [ψ(ξ[t+1] )] := E t
1 X i ψ(ξ[t] , ξt+1 ) 2ℓ 1 2ℓ−1 1 2ℓ−1
i=1 ℓ−1 2X i=1 ℓ−1 2X
∀ℓ ∈ N0 ,
2i−1 ψ(ξ[t] , ξt+1 )
∀ℓ ∈ N,
2i ψ(ξ[t] , ξt+1 )
∀ℓ ∈ N,
i=1
i , i ∈ [2ℓ ], are drawn independently from the marginal distribution of ξ (if t = 0) or where the samples ξt+1 1 the conditional distribution of ξt+1 given ξ[t] (if t > 0). By construction, we thus have
b ℓ [ψ(ξ[t+1] )] = 1 E b ℓ,o [ψ(ξ[t+1] )] + 1 E b ℓ,e [ψ(ξ[t+1] )]. E t 2 t 2 t 15
(17)
Definition 3.17 (MLMC Estimator). For any t ∈ [T − 1], select a rate parameter rt ∈ (0, 1) and truncation point Mt ∈ N, and let the random variable λt follow the truncated geometric distribution Geo(rt |Mt ). Assume that λt is independent of all other random objects, and set qt (ℓ) := P(λt = ℓ) for ℓ = 0, . . . , Mt b t (x), t ∈ [T ], using the following backward and t ∈ [T − 1]. We then construct a sequence of estimators H b recursion. We initialize the recursion by setting HT (x) := fT (ξT , x). Next, for any t ∈ [T − 1] we set 1 1 b λt ,o 1 b λt ,e λt b b Ht (x) := ht (x) − ht (x) − ht (x) , (18) qt (λt ) 2 2 where b bℓ H b t+1 (x) hℓt (x) := ft ξt , E ∀ℓ = 0, . . . , Mt t as well as ℓ,o b b b ℓ,e H b ℓ,o H b t+1 (x) b t+1 (x) := := hℓ,e (x) f ξ , E and h (x) f ξ , E t t t t t t t t
∀ℓ ∈ [Mt ].
b 0,o b b b We also set b h0,e t (x) := ht (x) := 0. The MLMC estimator for F (x) is then defined as F (x) := E[H1 (x)]. The MLMC estimator is built on a forest of n1 scenario trees that display a random branching structure with 2λt branches per stage-t node, where λt is sampled from a truncated geometric distribution. For this reason we will sometimes refer to λt as the log-branching factor of stage t. If Mt = ∞ for all t ∈ [T − 1] and the integrands ft map to R, then our MLMC estimator reduces to the untruncated MLMC estimator proposed in [61]. If additionally T = 2, then it further reduces to the estimator discussed in [66, § 3.2]. Working with random (log-)branching factors offers significant advantages. We will show that the bias of the MLMC estimator is of the same order as that of the SAA estimator with a large fixed branching factor 2Mt at each stage t ∈ [T −1]. Since the MLMC estimator is constructed on a scenario tree with a random branching factor 2λt that never exceeds but is typically much smaller than 2Mt at every stage, it achieves this bias using substantially fewer samples in total. As λt concentrates near zero, the branching factor 2λt often equals 1, so the MLMC estimator requires only a single sample in most branches. Consequently, the scenario tree grows more slowly with the time horizon T , making the evaluation of the MLMC estimator significantly more efficient than that of the SAA estimator with comparable bias. We now prove that the MLMC estimator requires fewer samples than the SAA estimator to achieve a mean squared error no greater than ϵ2 . To do so, we derive upper bounds on the variance, bias, and expected sampling cost of the MLMC estimator in Lemmas 3.19, 3.20, and 3.22, respectively. These bounds depend on the estimator’s hyperparameters, that is, the number n1 of scenario trees as well as the truncation points Mt and the rates rt that define the truncated geometric distributions of the log-branching factors λt , t ∈ [T − 1]. We then tune √ these hyperparameters to ensure that the bias and the variance of the MLMC estimator are bounded by ϵ/ 2. This calibration guarantees an overall mean squared error of at most ϵ2 , enabling a fair comparison of the expected sampling costs of the optimally tuned SAA and MLMC estimators. b t (x) conditional on time t − 1 information by From this point onward, we denote the p-th moment of H b t (x)∥p , µpt (x) := Et−1 ∥H (19) p where p ≥ 1. To analyze the MLMC estimator, we also define the conditional moment bound of order p as h i p µpT := sup ess sup ET −1 fT (ξT , x) p . (20) x∈Rd
′
By using Jensen’s inequality, one readily verifies that if µpT < ∞, then µpT < ∞ for every p′ ∈ [1, p]. The following assumption is crucial for our subsequent results. 16
Assumption 3.18. The conditional moment bound µpT of the integrand fT (ξT , x) is finite up to order p = 2 if Assumption 3.3 does not hold and up to order p = 2T if Assumption 3.3 holds. t
The following lemma establishes bounds on µ2t (x) and µ2t (x) for all t ∈ [T ], corresponding to the cases of nonsmooth and smooth integrands, respectively. In particular, for t = 1, these bounds can be used to derive upper bounds on the unconditional variance of the MLMC estimator in both settings. For any t ∈ [T ] Qt := and s ∈ [t + 1] we henceforth set L[s:t] i=s Li if s < t + 1; := 1 if s = t + 1; see also Lemma C.3 in the appendix. In particular, note that L[1:t] coincides with the Lipschitz constant L[t] defined before Lemma 3.6. Lemma 3.19. Suppose that Assumptions 3.1, 3.2 and 3.18 hold. Then, for all x ∈ X and t ∈ [T ], the b t (x) satisfy conditional moments of H Ms TY −1 X 1 2 if Assumption 3.3 does not hold, µt (x) ≤ Ct 2ℓ qs (ℓ) s=t ℓ=0 Ms TY −1 X 1 t µ2t (x) ≤ Dt if Assumption 3.3 holds, 2s ℓ q (ℓ)2s −1 2 s s=t ℓ=0
where Ct := µ2T (2B2 )T −t L2[t:T −1]
and
T Dt := µ2T
TY −1 s=t
3Ss 2
2s
s d2s −1 B2s+1
,
and where B2 and B2s+1 are universal constants defined in Lemma A.4 in the appendix. bH b 1 (x)]. Indeed, Lemma 3.19 allows us to bound the variance of the MLMC estimator Fb(x) = E[ b 1 (x), its variance admits the bound as Fb(x) constitutes an average of n1 independent copies of H V(Fb(x)) =
2 1 b 1 (x)) ≤ µ1 (x) , V(H n1 n1
(21)
while Lemma 3.19 provides a bound on µ21 (x). In fact, Lemma 3.19 constructs bounds on all conditional t moments µ2t (x) (for nonsmooth integrands that may violate Assumption 3.3) and µ2t (x) (for smooth integrands that satisfy Assumption 3.3) by backward recursion on t. Note that if the integrands are nonsmooth, then the upper bound on µ21 (x) is finite only if fT (ξT , x) has a finite second moment. If the integrands are smooth, on the other hand, then the upper bound on µ21 (x) is finite only if fT (ξT , x) has a finite moment of order 2T . The requirement that moments of order higher than two of fT (ξT , x) be finite arises as an artifact of the bound’s construction, which involves iterative application of the inequality (13) for smooth functions. This inequality is conservative because it neglects the Lipschitz continuity of fT (ξT , x) with respect to x. Indeed, the Lipschitz continuity ensures that the global growth of the left hand side of (13) is at most linear. In spite of this, we do not attempt to improve the conditional moment bounds for smooth integrands in Lemma 3.19. Deriving such bounds is cumbersome and unlikely to improve our scenario complexity results below. Moreover, to establish high-probability scenario complexity guarantees for the MLMC estimator, we will anyway assume that fT (ξT , x) is sub-Gaussian, which implies that all moments of fT (ξT , x) are finite. Lemma 3.19 readily extends to the untruncated MLMC estimator. In addition, if Assumption 3.3 holds bt while dt = 1 and Mt = ∞ for all t ∈ [T − 1], then our upper bound on the 2t -th conditional moment of H t coincides with the 2 -th unconditional moment of the untruncated MLMC estimator analyzed in [61]. The next result uses Lemma 3.19 to derive bounds on the bias of the MLMC estimator. Lemma 3.20. If Assumptions 3.1, 3.2 and 3.18 hold, then the bias of the MLMC estimator satisfies 1 PT −1 L[t] E[µ2t+1 (x)] 2 if Assumption 3.3 does not hold, 2Mt /2 E[Fb(x)] − F (x) ≤ Pt=1 T −1 L[t−1] St E[µ2t+1 (x)] if Assumption 3.3 holds. t=1 2Mt +1 17
Lemma 3.20 shows that the bias of the MLMC estimator decays faster with the truncation points if the integrands are smooth (i.e., if Assumption 3.3 holds) and if E[µ2t+1 (x)] is finite. Lemma C.6 in the appendix b t (x) conditional establishes a stronger result than Lemma 3.20 in that it provides a bound on the bias of H b b on time t − 1 information for every t ∈ [T ]. As E[F (x)] = E[H1 (x)], Lemma 3.20 follows immediately from Lemma C.6 for t = 1, and therefore its proof is omitted for brevity. Lemmas 3.20 and 3.6 reveal that the bias of the MLMC estimator with random branching factors 2λt for t ∈ [T − 1] is of the same order of magnitude as that of the SAA estimator with deterministic branching factors nt+1 = 2Mt for t ∈ [T −1]. In addition, Lemma 3.20 generalizes [41, Proposition 4.1], which focuses on conditional stochastic programs with T = 2, and [61, Theorem 2.2], which focuses on untruncated and thus unbiased MLMC estimators with Mt = ∞ for all t ∈ [T − 1]. Note that when T = 2, the bias of the SAA and MLMC estimators can be analyzed in the same manner. However, when T > 2, the bias of the MLMC estimator requires a more involved analysis due to its recursive construction, and it also involves a b t (x). bound on the variance of H Remark 3.21. The bias bounds stated in Lemma 3.20 are derived by first expressing the bias of Fb(x) as a telescoping sum. Then, Jensen’s inequality is applied to bound the absolute value of this sum by the sum of the absolute values of its terms. A tighter bias bound can be obtained using an inductive proof technique that, at each step, leverages a bias-variance decomposition. If the integrands are only known to be Lipschitz continuous (i.e., if Assumption 3.3 does not hold), then this refined technique yields the bias bound E[Fb(x)] − F (x) ≤
TX −1 L2 E[µ2 (x)] 1/2 t+1 [t] t=1
2Mt
. 1
This bound on the bias of Fb(x) matches the 2-norm of the vector with entries L[t] E[µ2t+1 (x)] 2 /2Mt /2 , t ∈ [T − 1]. Thus, it strengthens the corresponding bias bound of Lemma 3.20, which matches the 1-norm of the same vector. In addition, if the integrands are known to be smooth (i.e., if Assumption 3.3 holds), it is still possible to recursively construct a bias bound that improves upon the one given in Lemma 3.20. However, this refined bound does not admit a simple closed-form expression. As the recursively constructed bounds (for nonsmooth as well as smooth integrands) exhibit the same asymptotic dependence on the truncation points Mt , t ∈ [T − 1], we choose to work with the simpler and more interpretable bounds of Lemma 3.20. The next lemma provides an explicit formula for the expected sampling cost of the MLMC estimator. Lemma 3.22. The expected number of scenarios needed to construct the MLMC estimator amounts to Mt TY −1 X E C Fb(x) = n1 qt (ℓ)2ℓ . t=1 ℓ=0
Lemmas 3.19, 3.20 and 3.22 provide bounds on the variance, the bias and the expected computational cost of the MLMC estimator, respectively. They hold for any choices of the estimator’s hyperparameters, that is, the number n1 of scenario trees as well as the rate parameter rt and truncation point Mt of the geometric distribution governing the log-branching factor λt for every t ∈ [T − 1]. These hyperparameters should be tuned to guarantee a mean squared error of at most ϵ2 at minimal computational cost. To achieve this, we select sufficiently large truncation points Mt , t ∈ [T − 1], to ensure that the bias bound of Lemma 3.20 is √ at most ϵ/ 2. We also select a sufficiently large number n1 of scenario trees to ensure that the variance is at most ϵ2 /2. In view of (21), this is achieved by setting n1 ≈ 2µ21 (x)/ϵ2 , which implies that the expected sampling cost satisfies T −1 Mt 2µ2 (x) Y X C(Fb(x)) ≈ 12 qt (ℓ)2ℓ ; ϵ t=1 ℓ=0
18
see Lemma 3.22. Recall that we only have access to upper bounds on µ21 (x) that depend indirectly on the rate parameters rt through the probabilities qt (ℓ). We thus construct approximate upper bounds on C(Fb(x)) uniformly across all x ∈ X for both nonsmooth and smooth integrands by replacing µ21 (x) on the right hand side of the above expression with the respective bounds given in Lemma 3.19. We denote the resulting upper bounds by cns (r) and cs (r) for nonsmooth and smooth integrands, respectively, where r = (r1 , . . . , rT −1 ). Ideally, we should choose r ∈ (0, 1)T −1 so as to minimize these upper bounds. However, the exact minimizers are difficult to obtain in closed form. The next assumption provides near-optimal choices for r. Assumption 3.23 (Rate Parameters). For any t ∈ [T − 1] we set ( 2−1 if Assumption 3.3 does not hold, rt = −t −1−2 1−2 if Assumption 3.3 holds. Remark 3.24. Lemma C.7 in the appendix shows that if the integrands are nonsmooth (Assumption 3.3 does not hold), then the rate parameters of Assumption 3.23 exactly minimize the approximate upper bound cns (r) on the expected sampling cost over all r ∈ (0, 1). If Mt = ∞ for all t ∈ [T − 1], however, then the choice rt = 2−1 leads to infinite variance and cost bounds. In fact, no value of rt can simultaneously ensure finite second moment and finite expected sampling cost; see also [61]. These findings underscore the merits of truncated MLMC estimators. If the integrands are smooth (Assumption 3.3 holds), then the rate parameters of Assumption 3.23 exactly minimize the approximate upper bound cs (r) only for T = 2. For T > 2, the rate parameters in Assumption 3.23 jointly control both the geometric growth of the second moment and the expected sampling cost of Fb(x). Specifically, Lemma 3.19 implies that µ21 is bounded in Mt whenever t −1 rt < 1 − 2−1−(2 −1) , while Lemma 3.22 implies that the expected sampling cost is bounded in Mt whenever rt > 2−1 for all t ∈ [T − 1]. These bounds remain valid in the special case when Mt = ∞ for −t all t ∈ [T −1] addressed in [61]. The choice rt = 1−2−1−2 for all t ∈ [T −1] satisfies both of these strict inequalities and provides a balanced trade-off between the geometric growth of the level-ℓ contributions to the variance and the expected sampling cost. Hence, this choice is natural. The MLMC estimator is fully determined by the number n1 of scenario trees as well as the rates rt and the truncation points Mt that define the truncated geometric distributions of the log-branching factors λt for all t ∈ [T − 1]. The next theorem specifies a particular choice of these hyperparameters that guarantees the mean squared error of the MLMC estimator to be bounded by ϵ2 . It also characterizes the estimator’s mean squared error-based scenario complexity, which we define as the expected value of its sampling cost. Theorem 3.25 (Mean Squared Error-Based Scenario Complexity of the MLMC Estimator). Suppose that Assumptions 3.1, 3.2, 3.18 and 3.23 hold and that ϵ > 0 is a given tolerance. Set n1 = ⌈supx∈X 2µ21 (x)/ϵ2 ⌉, and define the truncation points Mt , t ∈ [T − 1], through the backward recursion √ Q √ T −1 √ √ Mt = 2 log2 2L[t] Ct+1 s=t+1 Zs Ms + 1 (T − 1)/ϵ if Assumption 3.3 does not hold and 1t 2s −1 √ − Mss+1 QT −1 2 2 1−2 2−t 2t Z Mt = log2 2L[t−1] St Dt+1 (T − 1)/(2ϵ) s=t+1 s − 1 1−2
2s
if Assumption 3.3 holds. Here, Ct and Dt are defined as in Lemma 3.19, and Zt := (1 − (1 − rt )Mt +1 )/rt . Then, the MLMC estimator achieves a root mean squared error of at most ϵ uniformly across all x ∈ X . Thus, the mean squared error-based scenario complexity of the MLMC estimator amounts at most to O(log(ϵ−1 )2(T −1) ϵ−2 ) if Assumption 3.3 does not hold and to O(ϵ−2 ) if Assumption 3.3 holds. 19
Theorems 3.9 and 3.25 show that the MLMC estimator substantially outperforms the SAA estimator in terms of mean squared error-based scenario complexity. Specifically, for nonsmooth integrands, the MLMC estimator achieves a scenario complexity of O(log(ϵ−1 )2(T −1) ϵ−2 ), representing a significant improvement over the O(ϵ−2T ) complexity of the SAA estimator. Similarly, for smooth integrands, the MLMC estimator improves the scenario complexity from O(ϵ−(T +1) ) for SAA to O(ϵ−2 ). The MLMC estimator of Definib ℓ,o tion 3.17, which discards all terms b hℓt , b hℓ,e t and ht beyond the truncation point Mt for each t ∈ [T − 1], also outperforms the untruncated MLMC estimator proposed in [61]. While the untruncated estimator achieves a scenario complexity of O(ϵ−2 ) for smooth integrands, it suffers from infinite scenario complexity for nonsmooth integrands. In contrast, our MLMC estimator avoids this divergence by using truncation. Note that while the scenario complexity of our MLMC estimator scales gracefully with 1/ϵ, it exhibits exponential dependence on T when the integrands are nonsmooth. As a result, the expected computational cost grows rapidly with T . This scaling behavior reflects the inherent difficulty of evaluating nested compositions of conditional expectations and nonsmooth functions. Nevertheless, for applications characterized by small values of T , as discussed in Section 2, the MLMC estimator can be computed to a moderate relative accuracy (e.g., 10% or 1%) at a practical computational cost. We emphasize, however, that Theorem 3.25 provides the first finite scenario complexity guarantee for problem (MCCE) with nonsmooth integrands that avoids exponential terms of the order ϵ−T . As these scenario complexity guarantees pertain only to expected errors, it is expedient to complement them with bounds that hold with high probability. As in Section 3.1, such bounds are available for estimators with finite second moments under a sub-Gaussianity assumption. b 1 (x)). The estimator H b 1 (x), which was introduced in DefiniAssumption 3.26 (Sub-Gaussianity of H 2 tion 3.17, is sub-Gaussian with variance proxy ζ > 0 for any fixed x ∈ X . Remark 3.27 below outlines easily verifiable conditions under which Assumption 3.26 holds. Remark 3.27. Assumption 3.26 holds under the same conditions that imply Assumption 3.11 in Remark 3.12. b t (x) can be shown by backward recursion on t ∈ [T −1], using the assumed Indeed, the sub-Gaussianity of H sub-Gaussianity of ξt together with Lemma C.5, which ensures that sub-Gaussianity is preserved under finite sums, concatenations, Lipschitz continuous transformations and mixtures. Using Assumption 3.26, we now demonstrate that if the rates rt and truncation points Mt for all t ∈ [T −1] are chosen as in Theorem 3.25, then the worst-case absolute error of the MLMC estimator is bounded by ϵ with high probability, where the failure probability tends to 0 as n1 increases. In the remainder of this section, we use L′ ∈ R to denote the Lipschitz constant of the MLMC estimator Fb(x) with respect to x, which can be shown to exist by using similar arguments as in Lemma C.3. Details are omitted for brevity. Lemma 3.28. Suppose that Assumptions 3.1, 3.2, 3.23 and 3.26 hold and that ϵ > 0. Assume also that the truncation points Mt , t ∈ [T − 1] satisfy Q √ T −1 √ √ Mt = 2 log2 4L[t] Ct+1 s=t+1 Zs Ms + 1 (T − 1)/ϵ if Assumption 3.3 does not hold and 1t 2s −1 − Mss+1 QT −1 −t 2 2 1−2 2 2t Mt = log2 2L[t−1] St Dt+1 (T − 1)/ϵ s=t+1 Zs − 1 1−2
2s
if Assumption 3.3 holds. Then, we have ′ d n1 ϵ 2 b P sup |F (x) − F (x)| > ϵ ≤ 2 8L DX /ϵ + 1 exp − . 32ζ 2 x∈X 20
b 1 (x) (see Assumption 3.26), while d and DX are the dimension Recall that ζ 2 is the variance proxy of H and the diameter of X , respectively. The proof of Lemma 3.28 parallels that of Lemma 3.13 and is thus omitted. Lemma 3.28 readily implies that the minimizers of problem (1), which optimizes the MLMC estimator over X , are near-optimal in (MCCO) with high probability. Unfortunately, these minimizers are often difficult to compute even if F (x) is a convex function. Remark 3.29. Recall that F (x) as well as the SAA estimator from Section 3.1 are convex if the integrand ft (ξt , xt ) is convex in xt for every t ∈ [T ] and non-decreasing in xt for every t ∈ [T − 1]. Perhaps surprisingly, however, these conditions fail to guarantee that the MLMC estimator is convex. To see this, assume for example that T = 2, f1 (ξ1 , x1 ) = exp(x1 ) and f2 (ξ2 , x2 ) = ξ2 + x2 . As ft (ξt , xt ) is convex and non-decreasing in xt for every t ∈ [T ], we conclude that F (x) and its SAA estimator are convex. Next, we construct the MLMC estimator with n1 = 1 scenario trees in the event λ1 = 1. This estimator requires a sample ξ11 from the marginal distribution of ξ1 and 2λ1 = 2 samples ξ21 and ξ22 from the conditional distribution of ξ2 given ξ1 = ξ11 . By Definition 3.17, the estimator can be represented as 1 1 i 1 h exp (ξ21 + ξ22 )/2 + x − exp ξ21 + x − exp ξ22 + x q1 (λ1 ) 2 2 h 1 1 i exp(x) = exp (ξ21 + ξ22 )/2 − exp ξ21 − exp ξ22 , q1 (1) 2 2
Fb(x) =
where the second equality holds because we condition on the event λ1 = 1. Assuming that ξ21 ̸= ξ22 , Jensen’s inequality implies that the last term in square brackets is negative, which in turn implies that Fb(x) is nonconvex. Hence, minimizing the MLMC estimator over X may, with positive probability, lead to an intractable nonconvex optimization problem. We emphasize, however, that the MLMC estimator constitutes an unbiased estimator for the SAA estimator with branching factor 2M1 , which is convex in this example. We now demonstrate that the minimizers of the MLMC problem (1) are, regardless of whether they can be computed efficiently, near-optimal in (MCCO) with high probability. Corollary 3.30. Assume that all conditions of Lemma 3.28 hold for some ϵ > 0, and let x∗ and x b∗ be minimizers of (MCCO) and the MLMC problem (1), respectively. Then, we have ′ d n1 ϵ2 ∗ ∗ P F (b x ) − F (x ) > ϵ ≤ 4 8L DX /ϵ + 1 exp − . 128ζ 2 The proof of Corollary 3.30 parallels that of Corollary 3.15 and is thus omitted. As the computational cost of the MLMC estimator is random, we generalize our definition of the high-probability scenario complexity as the expected number of scenarios used by a method that is guaranteed to find an ϵ-optimal solution of problem (MCCO) with probability 1 − β. The next result follows immediately from Corollary 3.30. Theorem 3.31 (High-Probability Scenario Complexity of the MLMC Estimator). Suppose that all conditions of Lemma 3.28 hold for some ϵ > 0 and that ′ 128ζ 2 n1 = d log 8L D /ϵ + 1 + log (4/β) X ϵ2 for some β ∈ (0, 1). Then any solution of the MLMC problem (1) is an ϵ-optimal solution of (MCCO) with probability at least 1 − β. Thus, the high-probability scenario complexity of the MLMC estimator equals O(log(ϵ−1 )T ϵ−2 ) if Assumption 3.3 does not hold and O(log(ϵ−1 )ϵ−2 ) if Assumption 3.3 holds.
21
Theorems 3.16 and 3.31 show that the MLMC estimator significantly outperforms the SAA estimator in terms of high-probability scenario complexity. Indeed, it improves the high-probability scenario complexity from O(log(ϵ−1 )ϵ−2T ) to O(log(ϵ−1 )T ϵ−2 ) for nonsmooth integrands and from O(log(ϵ−1 )ϵ−(T +1) ) to O(log(ϵ−1 )ϵ−2 ) for smooth integrands. In contrast, the untruncated MLMC estimator studied in [61] b t (x) scales lacks a finite Lipschitz constant. Indeed, as is evident from (18), the Lipschitz modulus of H with 1/q(λt ). As the random variable λt has unbounded support when Mt = ∞, 1/q(λt ) can be arbitrarily b t (x) admits no Lipschitz constant that applies uniformly across all realizations of λt . Our large. Hence, H bound on the high-probability scenario complexity of the untruncated MLMC estimator thus diverges. Note that the high-probability guarantees of Theorem 3.31 for smooth integrands require a larger number n1 of scenario trees than the mean squared error-based scenario complexity results in Theorem 3.25. Specifically, the high-probability scenario complexity incurs an additional multiplicative factor of log(ϵ−1 ). Despite this increase, the ϵ-scaling of the MLMC estimator remains better than that of the SAA estimator. For non-smooth integrands, the high-probability scenario complexity O(log(ϵ−1 )T ϵ−2 ) even grows slower with ϵ−1 than the mean squared error-based sample complexity O(log(ϵ−1 )2(T −1) ϵ−2 ) from Theb 1 (x) in orem 3.25. This phenomenon arises because of the sub-Gaussian tail assumption on the estimator H Theorem 3.31, which yields tighter control over the worst-case behavior than variance-based bounds. Finally, we emphasize that while the scenario complexity of the truncated MLMC estimator grows gracefully with ϵ−1 for (moderate) fixed values of T , it usually grows exponentially with T . Such exponential growth appears to be unavoidable even for smooth integrands, in which case the scenario complexity can be shown to include a constant term that is exponential in T . In addition, the untruncated MLMC estimator from [61] as well as the SAA estimator from Section 3.1 also display an exponential cost growth in T .
4
Gradient Methods for MCCO
From Section 3 we know that the minimizers of problem (1), equipped with the SAA or the MLMC estimator, are near-optimal in (MCCO) with high probability. Despite their favorable statistical properties, computing these minimizers can be challenging in the absence of restrictive convexity and monotonicity conditions. This motivates us in this section to develop efficient estimators for ∇F (x) ∈ Rd that can be employed within first-order optimization methods for solving problem (MCCO) to stationarity. We aim to recursively construct MLMC gradient estimators using similar ideas as in Section 3.2. These estimators are again based on a forest of scenario trees with a random branching structure, and the underlying log-branching factors are modeled by truncated geometric distributions. The resulting estimators are biased, but the bias can be reduced by increasing the truncation points. In the special case where the truncation points are set to infinity, we obtain untruncated, and therefore unbiased, estimators. Both the truncated and untruncated gradient estimators proposed here are novel for problems with more than two stages (T > 2). Truncated as well as untruncated MLMC gradient estimators for two-stage problems are discussed in [29, 41]. As will become clear below, the generalization to T ≥ 2 poses significant technical challenges. Throughout, we work under Assumptions 3.2 and 3.3, which guarantee the existence and well-definedness of the gradients. We first show that the gradient ∇F (x) can be constructed via the backward recursion T Gt (x) := Et−1 Gt+1 (x)∇xt ft (ξt , Ft+1 (x)) ∈ Rd×dt−1 ∀t ∈ [T − 1] (22) T (x) represents the nest of the last initialized by GT (x) := ET −1 [∇x fT (ξT , x)] ∈ Rd×dT −1 , where Ft+1 T − t conditional expectations and nonlinear integrands in (MCCO). Hence, we have F (x) = F1T (x); see also Definition C.1 in the appendix. All gradients in (22) exist by Assumption 3.3. Moreover, since ∇xt ft (ξt , xt ) is uniformly bounded under Assumption 3.2, the same holds for Gt (x). Consequently, the conditional expectation in (22) is well-defined. We can now prove that ∇F (x) = G1 (x).
Lemma 4.1. If Assumptions 3.1, 3.2 and 3.3 hold, then ∇F (x) = G1 (x) for all x ∈ Rd . 22
Lemma 4.1 shows that the gradient ∇F (x) admits a recursive representation analogous to that of F (x) itself. This observation allows us to construct gradient estimators using methods similar to those in Section 3. In particular, following the analogy with Definition 3.17, we can define an MLMC estimator for ∇F (x). Definition 4.2 (MLMC Gradient Estimator). For any t ∈ [T − 1], select a rate parameter rt ∈ (0, 1) and truncation point Mt ∈ N, and let λt be a random variable that follows the truncated geometric distribution Geo(rt |Mt ). Assume that λt is independent of all other random objects, and set qt (ℓ) := P(λt = ℓ) b t (x), t ∈ [T ], using the following for ℓ = 0, . . . , Mt and t ∈ [T − 1]. We construct gradient estimators G b backward recursion. We initialize the recursion by GT (x) := ∇x fT (ξT , x). Next, for any t ∈ [T − 1], we set 1 1 1 λ ,e λ ,o λ t t t b t (x) := G gb (x) − gbt (x) − gbt (x) , (23) qt (λt ) t 2 2 where as well as
b ℓ [H b ℓ [G b t+1 (x)] b t+1 (x)] ∇xt ft ξt , E gbtℓ (x) := E t t
∀ℓ = 0, . . . , Mt
b ℓ,e [G b ℓ,e [H b t+1 (x)] ∇xt ft ξt , E b t+1 (x)] ∀ℓ ∈ [Mt ], gbtℓ,e (x) :=E t t b ℓ,o [G b ℓ,o [H b t+1 (x)] ∇xt ft ξt , E b t+1 (x)] ∀ℓ ∈ [Mt ]. gbtℓ,o (x) :=E t t
bG b b 1 (x)]. := E[ We also set gbt0,e (x) := gbt0,o (x) := 0. The MLMC estimator for ∇F (x) is then defined as G(x) A special case arises when the truncation points in Definition 4.2 are set to Mt = ∞ for all t ∈ [T −1]. In b t (x) this case, we refer to the resulting estimator as the untruncated MLMC gradient estimator. Recall that H was introduced in Definition 3.17 and can be viewed as an estimator for the conditional risk mapping FtT (x) b t (x) constitutes an estimator for Gt (x). Crucially, the construction above from Definition C.1. Similarly, G b t (x) depends on the function reveals an important structural feature, namely that the gradient estimator G b value estimator Ht+1 (x) as an argument at each stage, and both sequences of estimators must therefore be maintained and propagated simultaneously. This coupling is a novel structural feature of our construction, b t (x) and H b t (x) are built on the same forest of n1 absent from classical MLMC gradient estimators. Both G λ t scenario trees displaying a random branching structure with 2 branches per stage-t node, where λt is sampled from a truncated geometric distribution. Note that the forest of trees is defined implicitly through b and E b λt [·] for t ∈ [T − 1], which are used both the nesting of the empirical conditional expectations E[·] t b t (x) and H b t (x), for the construction of the function value and gradient estimators. This implies that G t ∈ [T − 1], are constructed from the same samples and are therefore correlated. Reusing the same samples across both estimators reduces the sampling cost, but it complicates the analysis of the trade-off between variance and computational cost. Our analysis uses the following Hölder condition. Assumption 4.3 (Hölder Continuous Hessians). The Hessian matrix ∇2xt ft,i (ξt , xt ) is ρt -Hölder continuous in xt across all ξt . That is, for all t ∈ [T − 1], there exist Rt > 0 and ρt ∈ (1 − 21−t , 1] such that ∥∇2xt ft,i (ξt , xt ) − ∇2xt ft,i (ξt , yt )∥2 ≤ Rt ∥xt − yt ∥ρ2t
∀xt , yt ∈ Rdt , ∀ξt ∈ Rmt , ∀i ∈ [dt−1 ].
Note that when ρt = 1, Assumption 4.3 requires the Hessian matrices of the components of ft (ξt , xt ) to be Lipschitz continuous in xt . Assumption 4.3 has been used in [29] to control the variance and the sampling costs of MLMC gradient estimators for conditional stochastic optimization problems, which arise as special cases of (MCCO) with T = 2 nests. It implies via the fundamental theorem of calculus that ∇xt ft,i (ξt , xt ) − ∇xt ft,i (ξt , yt ) − ∇2xt ft,i (ξt , yt )(xt − yt ) 2 ≤ for all xt , yt ∈ Rdt and ξt ∈ Rmt . 23
Rt ∥xt − yt ∥ρ2t +1 ρt + 1
(24)
In order to derive bounds on the mean squared error-based scenario complexity of the MLMC gradient estimator, we first analyze its variance, bias, and sampling cost. In analogy with (19), for any p ≥ 1, we b t (x) given the information available up to time t − 1 by define the p-th conditional moment of G h i b t (x)∥p . νtp (x) := Et−1 ∥G (25) 2 b t (x)∥2 stands for its Frobenius norm. Note also that b t (x) ∈ Rd×dt−1 is a matrix, and thus ∥G Recall that G b the p-th conditional moment of Ht (x) in (19) was defined with respect to the vector p-norm. To facilitate the analysis of the MLMC gradient estimator, we define the conditional moment bound of order p as i h p ν pT := sup ess sup ET −1 ∇x fT (ξT , x) p . x∈Rd
The following assumption is crucial for our subsequent results. Assumption 4.4. The conditional moment bound µpT of the integrand fT (ξT , x) is finite up to order p = 2T (ρT −1 +1), and the conditional moment bound ν pT of the gradient ∇fT (ξT , x) is finite up to order p = 2T . Note that Assumption 4.4 implies Assumption 3.18 because ρT −1 > 0. The following lemma establishes t a bound on νt2 (x) for all t ∈ [T ]. In particular, for t = 1, this bound can be used to derive an upper bound on the unconditional variance of the MLMC gradient estimator. Lemma 4.5. Suppose that Assumptions 3.1, 3.2, 3.3, 4.3 and 4.4 hold. Then, for every t ∈ [T − 1] there exists Et that depends only on T , the moment bounds from Assumption 4.4, the dimensions ds and smoothness constants Ss for s = t + 1, . . . , T − 1, and the Hölder constants Rs and ρs for s = t, . . . , T such that t νt2 (x) ≤ Et
TY −1 s=t
max
X Ms
Ms X
1 , s s s−1 s 2(2 )(ρs +1)ℓ qs (ℓ)2 −1 ℓ=0 22 (ρs−1 +1)ℓ qs (ℓ)2 (ρs−1 +1)−1 ℓ=0 1
for all x ∈ X and t ∈ [T ]. Here, we adopt the convention ρ0 = 0 for notational convenience. The exact closed-form expressions for the constants Et are cumbersome and thus omitted from the statement of Lemma 4.5. Instead, an explicit recursive definition is given in the proof. Lemma 4.5 allows us bG b b 1 (x)]. Indeed, as G(x) b to bound the variance of the MLMC gradient estimator G(x) = E[ constitutes an b 1 (x), its variance admits the bound average of n1 independent copies of G b V(G(x)) =
2 1 b 1 (x)) ≤ ν1 (x) , V(G n1 n1
(26)
while Lemma 4.5 provides a bound on ν12 (x). Note that this bound is finite only if ∇fT (ξT , x) admits finite moments of order 2T , and fT (ξT , x) admits finite moments of order 2T (ρT −1 + 1). We further emphasize that Lemma 4.5 also applies to the untruncated MLMC gradient estimator obtained by setting Mt = ∞ for all t ∈ [T − 1]. In this case, however, the resulting variance bound remains finite only if the rate parameters rt , t ∈ [T − 1], are sufficiently large to ensure convergence of all infinite series appearing in the bound. Lemma 4.5 is closely related to Lemma 3.19, as both establish bounds on conditional moments of MLMC estimators by means of backward recursions in time. The key differences lie in the underlying assumptions and in the complexity of the resulting recursions. For function-value estimators, Lemma 3.19 requires only Lipschitz continuity of the integrands, yielding relatively simple product-form bounds that admit closed-form expressions. If the integrands are also smooth, the bounds strengthen while retaining a simple closed form. By contrast, for gradient estimators, Lemma 4.5 requires not only smoothness of the integrands but also Hölder continuity of their Hessian matrices, which introduces the Hölder exponents ρt 24
into the analysis and makes the proof technically more demanding. Nevertheless, both lemmas serve the same purpose. They provide explicit upper bounds on conditional moments that can be used to control the growth of the estimator’s variance and to establish scenario complexity guarantees. When ρt = 1 for all t ∈ [T − 1] (i.e., when all Hessian matrices are Lipschitz continuous), the second term in the max operator of the moment bound in Lemma 4.5 dominates the first term for all s ∈ [T − 1]—recall that ρ0 = 0—which implies that the conditional moment bounds for the gradient estimator impose stricter requirements on the rate parameters rs governing qs (ℓ) than those for the function-value estimator in Lemma 3.19. Next, we analyze the bias of the MLMC gradient estimator using the moment bounds established in Lemmas 3.19 and 4.5. This analysis is cumbersome because gbtℓ and its variants based on only even- or odd-indexed samples, gbtℓ,e and gbtℓ,o , respectively (see Definition 4.2), exhibit a product structure. More b t+1 is propagated through the gradient mapping ∇xt ft (ξt , ·), and its magnitude is precisely, any bias in H b t+1 and H b t+1 interact multiplicatively, thus controlled by the smoothness of ft . Furthermore, the biases of G giving rise to cross terms that must be carefully controlled. This interplay significantly complicates the analysis compared with the case of the function value estimators studied in Section 3.2. Lemma 4.6. If Assumptions 3.1, 3.2, 3.3 and 4.4 hold, the bias of the MLMC gradient estimator satisfies 1 1 T −1 T −1 X X L[t+1:s−1] Ss E[(µ2s+1 (x))2 ] 2 E[µ2t+1 (x)] 2 1 2 b E[G(x)] − G1 (x) 2 ≤ L[t−1] St E[νt+1 (x)] 2 . + 2Ms +1 2Mt /2 t=1
s=t+1
As in the case of the MLMC function value estimator discussed in Section 3.2, Lemma 4.6 shows that the MLMC gradient estimator is unbiased if Mt = ∞ for all t ∈ [T − 1] provided that all expectations in the bias bound are finite. Hence, the untruncated MLMC gradient estimator is unbiased. We now quantify the expected sampling cost of the MLMC gradient estimator. Lemma 4.7. The expected number of scenarios needed to construct the MLMC gradient estimator equals Mt TY −1 X b E C G(x) = n1 qt (ℓ)2ℓ . t=1 ℓ=0
The proof of Lemma 4.7 widely parallels that of Lemma 3.22 and is thus omitted. As in Section 3.2, we now tune the hyperparameters n1 , rt and Mt of the MLMC estimator to guarantee a mean squared error of at most ϵ2 at minimal computational cost. To achieve this, we select sufficiently large truncation points Mt , √ t ∈ [T − 1], to ensure that the bias bound of Lemma 4.6 is at most ϵ/ 2. We also select a sufficiently large number n1 ≈ 2ν12 (x)/ϵ2 of scenario trees to ensure that the variance is at most ϵ2 /2; see (26). By Lemma 4.7, the expected sampling cost can thus be approximated by T −1 M
b C(G(x)) ≈
t 2ν12 (x) Y X qt (ℓ)2ℓ . ϵ2
t=1 ℓ=0
b We can construct an approximate upper bound on C(G(x)) uniformly across all x ∈ X by replacing ν12 (x) with the respective bound given in Lemma 4.5. The resulting upper bound ch (r) depends on the rate parameters r = (r1 , . . . , rT −1 ) of the truncated geometric distributions used to model the log-branching factors. Here, the subscript “h” indicates that the bound was obtained under the assumption of Hölder continuous Hessians. Ideally, r ∈ (0, 1)T −1 should be selected to minimize the given upper bound. As in Section 3.2, closed-form minimizers are elusive. The next assumption specifies a simple near-optimal choice for r. Assumption 4.8 (Rate Parameters). For any t ∈ [T − 1], we select rt satisfying 2t (ρ +1) 2t−1 (ρt +1) − 2t (ρ t−1 − −1 +1)−1 t 2 −1 t−1 , 1−2 . 2 < rt < min 1 − 2 25
A rate parameter rt that satisfies both strict inequalities is guaranteed to exist whenever Assumption 4.3 holds. Indeed, in this case the Hölder exponent satisfies ρt > 1 − 21−t , which guarantees that the first term in the minimum is strictly greater than 2−1 . The second term is strictly greater than 2−1 for any ρt ≥ 0. Assumption 4.8 ensures that the given rate parameters jointly control the geometric growth of the second b moment as well as the expected sampling cost of G(x). Specifically, Lemma 4.5 implies that ν12 remains bounded in Mt if rt obeys the given upper bound in Assumption 4.8, and Lemma 4.7 implies that the expected sampling cost is bounded in Mt if rt > 2−1 for all t ∈ [T − 1]. Note also that if Mt = ∞ for all t ∈ [T − 1], then Assumption 4.8 guarantees that the infinite sums in Lemmas 4.5 and 4.7 converge. Hence, the variance and the expected sampling cost of the untruncated MLMC gradient estimator are finite. The MLMC gradient estimator is fully characterized by the number n1 of scenario trees as well as the rates rt and truncation points Mt governing the geometric distributions of the log-branching factors λt for t ∈ [T − 1]. The following theorem presents a specific configuration of these hyperparameters that ensures the root mean squared error of the truncated MLMC gradient estimator to be bounded by ϵ. Theorem 4.9 (Mean Squared Error-Based Scenario Complexity of the MLMC Gradient Estimator). Let Assumptions 3.1, 3.2, 3.3, 4.3, 4.4 and 4.8 hold, and let ϵ > 0 be a given tolerance. Then, for every t ∈ [T − 1] there exists Wt that depends only on T , the moment bounds from Assumption 4.4, the dimensions ds and smoothness constants Ss for s = t + 1, . . . , T − 1, the Lipschitz constants Ls and the rate parameters rs for s ∈ [T − 1], and the Hölder constants Rs and ρs for s = t, . . . , T such that the following holds. If n1 = ⌈sup 2ν12 (x)/ϵ2 ⌉ x∈X
and
Mt = ⌈2 log2 (Wt /ϵ)⌉
∀t ∈ [T − 1],
then the MLMC gradient estimator achieves a root mean squared error of at most ϵ uniformly for all x ∈ X . Thus, its mean squared error-based scenario complexity is at most O(ϵ−2 ). The exact closed-form expressions for the constants Wt are cumbersome and thus omitted from the statement of Theorem 4.9. The explicit construction of these constants will become clear in the proof. Theorem 4.9 establishes that, if the integrands have Hölder continuous Hessians, then the MLMC gradient estimator achieves the optimal scenario complexity of O(ϵ−2 ). This may be surprising because the MLMC gradient estimator has a significantly more complex structure than the MLMC function value estimator, implying that it is more difficult to control its variance; see also the discussion after Lemma 4.5. We emphasize that Theorem 4.9 remains valid for the untruncated MLMC gradient estimator, which is obtained by setting Mt = ∞ for all t ∈ [T − 1]. This result is new for T > 2. When T = 2, Theorem 4.9 recovers the O(ϵ−2 ) scenario complexity for estimating the gradients of conditional stochastic optimization problems established in [29, Theorem 2]. The Hölder continuity of the Hessians is a crucial assumption for Theorem 4.9. Indeed, if the integrands are merely Lipschitz continuous and smooth, then no choice of rt , t ∈ [T − 1], guarantees that the variance and the expected sampling cost remain simultaneously bounded [41]. Hence, the scenario complexity of the untruncated MLMC gradient estimator diverges. Remark 4.10. The MLMC gradient estimator introduced in Definition 4.2 employs the same set of samples b t (x) and H b t (x). This shared use of samples reduces the expected sampling cost. Howto construct both G b t (x) and H b t (x), which complicates the analysis of the scenario ever, it also induces correlation between G complexity and necessitates the assumption that the integrands possess Hölder continuous Hessians. Alternatively, one could construct a simpler MLMC gradient estimator that employs independent sets of samples b t (x) and H b t (x). The independence between these components simplifies the to compute the components G moment bounds in Lemma 4.5 and, in fact, eliminates the need to establish moment bounds of orders p > 2. Consequently, it suffices to derive a second-order moment bound analogous to that obtained in Lemma 3.19 for nonsmooth integrands. However, this alternative estimator entails a higher expected sampling cost. One can show that its mean squared error-based scenario complexity amounts to O(log(ϵ−1 )2(T −1) ϵ−2 ), which 26
matches that of the function value estimator for nonsmooth integrands from Theorem 3.25. The proof follows similar arguments to those used in Theorem 4.9, albeit in a simpler form. Detailed derivations are therefore omitted for brevity. We emphasize that, unless the integrands have Hölder continuous Hessians, the truncation points of the simpler MLMC gradient estimator cannot be sent to infinity without causing either the variance or the expected sampling cost to diverge. Hence, in this simpler setting, it is impossible to construct an untruncated unbiased estimator without imposing additional regularity conditions. We can now show that the MLMC gradient estimator introduced in Definition 4.2 can be incorporated into an SGD algorithm (see Algorithm to efficiently compute an ϵ-stationary point of problem (MCCO). P1) T 2 := In the following, we use S[T ] t=1 L[1:t−1] St L[t+1:T ] as a measure for the smoothness of F (x), where L[s:t] for t ∈ [T ] and s ∈ [T + 1] is defined as before Lemma 3.19. For more details see Lemma C.3. 2 2 Theorem 4.11. Let Assumptions 3.1, 3.2, 3.3, 4.3, p 4.4 and 4.8 hold, and set ν 1 := supx∈X ν1 (x). Assume that Algorithm 1 uses a constant stepsize η = ν 1 n1 /K and uses the MLMC gradient estimator with n1 = O(1) and Mt , t ∈ [T − 1], defined as in Theorem 4.9. If Algorithm 1 is initialized at x1 ∈ X and runs over
i 2 ν 21 h K≥ 2 F (x1 ) − min F (x) + S[T ] x∈X n1 ϵ4 iterations, then it outputs an ϵ-stationary point x̄K with E[∥∇F (x̄K )∥22 ] ≤ ϵ2 . The last part of the proof of Theorem 4.9 shows that if n1 = O(1), then the expected sampling cost of constructing a single MLMC gradient estimator is independent of ϵ (i.e., O(1)). Theorem 4.11 thus implies that Algorithm 1 requires O(ϵ−4 ) scenarios to compute an ϵ-stationary point of (MCCO). Theorem 4.11 generalizes a known scenario complexity result for conditional stochastic optimization problems, which correspond to instances of (MCCO) with T = 2 stages [41, Corollary 4.1]. When employing untruncated MLMC gradient estimators, obtained by setting Mt = ∞ for all t ∈ [T − 1], a general complexity result for constrained nonconvex stochastic optimization with unbiased gradient oracles [22, Corollary 4] implies that an ϵ-stationary point of problem (MCCO) can still be computed in O(ϵ−4 ) iterations. This convergence rate is known to be optimal for unbiased gradient oracles in the absence of variance reduction [2]. Variancereduction techniques can further improve the optimal rate for unbiased gradient oracles to O(ϵ−3 ) [2]. Algorithm 1 Projected SGD Input: Number of iterations K, stepsize schedule {ηk }K k=1 , initial iterate x1 for k = 1, . . . , K do b k ) for the gradient ∇F (xk ) Construct an estimator G(x b k )), where ΠX is the Euclidean projection on X Update xk+1 = ΠX (xk − ηk G(x end for Output: x̄K chosen uniformly at random from {xk }K k=1
5
Numerical Experiments
We now assess the benefits of our MLMC function value and gradient estimators for solving MCCE and MCCO problems, and we experimentally validate our theoretical results. All experiments are conducted on a MacBook Pro equipped with an Apple M1 Max chip and 32 GB RAM, running macOS Sonoma 14.5. Code and data for reproducing the numerical experiments are available at https://github.com/busesen/ MCCO. 27
MCCE with a Known Ground-Truth Solution The first experiment is centered around a synthetic instance of (MCCE), originally introduced in [61, § 3.1], with T = 3 stages and integrands f1 (ξ1 , x1 ) = sin(ξ1 + x1 ), f2 (ξ2 , x2 ) = sin(ξ2 − x2 ) and f3 (ξ3 , x) = ξ3 . The disturbances follow a Markovian structure with ξ1 ∼ N (π/2, 1), ξ2 |ξ1 ∼ N (ξ1 , 1) and ξ3 |ξ2 ∼ N (ξ2 , 1). This problem has the closed-form solution F (x) = exp(−1/2) ≈ 0.6065, against which the accuracy of different estimators can be assessed. We compare two SAA estimators in the sense of Definition 3.4 with two MLMC estimators in the sense of Definition 3.17. The first SAA estimator (SAA1) employs a uniform branching structure with n1 = n2 = n3 , while the second SAA estimator (SAA2) emphasizes the first stage by setting n1 = n22 = n23 ; see [51, 61] for motivation and discussion. The first MLMC estimator, taken from [61], is untruncated, uses the recommended rate parameters r1 = 0.74 and r2 = 0.60, and sets M1 = M2 = ∞, resulting in an expected sampling cost of 4.6250 per each of the n1 underlying scenario trees. The second MLMC estimator is truncated, adopts the rates r1 = 1 − 2−3/2 and r2 = 1 − 2−5/4 proposed in Assumption 3.23 for smooth integrands, and sets M1 = 6 and M2 = 5, yielding an expected sampling cost of 4.7674 per scenario tree. The truncation points are chosen to satisfy M1 > M2 , in accordance with Lemma 3.28, while maintaining an expected sampling cost per tree comparable to that of the untruncated MLMC estimator. For each of the four estimators, we consider a sequence of instances indexed by n1 , yielding decreasing mean squared errors as n1 increases. For the SAA estimators, the associated scenario trees become increasingly bushy, and arguments analogous to those in Theorem 3.9 show that their mean squared errors converge to 0, implying asymptotic unbiasedness. By Lemma 3.20, the untruncated MLMC estimator is unbiased for every fixed n1 , and similar reasoning as in Theorem 3.25 establishes that its mean squared error, which coincides with its variance, also vanishes as n1 grows. In contrast, the truncated MLMC estimator exhibits a nonvanishing bias as n1 grows. Thus, increasing n1 serves only to reduce its variance. Figure 2 (a) displays the two MLMC estimators (lines) together with their 95% confidence intervals (shaded regions) as functions of n1 . Since any MLMC estimator Fb(x) is formed by averaging n1 independent tree-wise estimators with finite second moments (see Definition 3.17 and Lemma 3.19), the central limit theorem implies approximate normality, and we therefore report confidence bounds of the form √ Fb(x) ± 1.96 σ bn1 / n1 , where σ bn1 denotes the empirical standard deviation of the n1 tree-wise estimators. For n1 < 500 the variance becomes excessively large, and the plot is therefore restricted to n1 ≥ 500. Although unbiased, the untruncated MLMC estimator exhibits substantially larger variance across all values of n1 and is notably fragile, as illustrated by the discontinuity at n1 = 17,212. Such jumps are caused by the occasional sampling of large log-branching factors, which can occur because the underlying geometric distribution is not truncated.1 The figure shows the first of 10 independent simulation runs, with similarly pronounced discontinuities observed in 4 out of the 10 runs for the untruncated estimator. Overall, truncation yields markedly more stable MLMC estimators with reduced variance, and although it introduces bias, this bias remains negligible at least up to n1 = 105 . The two SAA estimators are known to converge slower to F (x) than the MLMC estimators and are thus not shown in Figure 2 (a); see, e.g., [61, § 3]. Figure 2 (b) illustrates the trade-off between mean squared error and sampling cost for the four estimators. As before, different instances of each estimator are parametrized by n1 . Since the truncated MLMC estimator is biased, its mean squared error is eventually dominated, for sufficiently large n1 , by the squared bias. This causes the error to saturate as n1 (and thus the total sampling cost) increases. By Lemma 3.20, the squared bias is bounded above by (2−7 + 2−6 )2 ≈ 5.4932 × 10−4 , implying that saturation occurs below a log-mean squared error of log10 (5.4932 × 10−4 ) ≈ −3.2602. It actually lies outside of the plotting range. The untruncated MLMC estimator is unbiased for any n1 , yet it attains a larger mean squared error for all shown computational budgets due to its substantially higher variance. Although for sufficiently large n1 this variance can be reduced below the squared bias of the truncated MLMC estimator, so that the 1
The discontinuity in Figure 2 (a) occurs because a large log-branching factor equal to ℓ = 10 was sampled in the 17,212-th tree, which incurs a sampling cost that scales with 2ℓ and also introduces a large importance weight 1/qt (ℓ) in (18).
28
untruncated estimator may eventually outperform its truncated counterpart with fixed truncation levels, the required sampling costs typically exceed practical computational limits. Moreover, with sufficiently large computational budgets, one could jointly increase n1 as well as the truncation points, thereby reducing both variance and bias of the truncated estimator and preserving its advantage. The mean squared errors of the two MLMC estimators decay at comparable empirical rates of −0.8707 and −0.7871 for the truncated and untruncated variants, respectively, in agreement with the theoretical n−1 1 rate for smooth integrands (see Theorem 3.25 and [61, Corollary 2.3]). In contrast, the SAA estimators perform substantially worse, with −1/3 decay rates of −0.3352 for SAA1 and −0.4936 for SAA2, consistent with the theoretical rates n1 and −1/2 n1 , respectively (see Theorem 3.9 and [51, Appendix G]). Taken together, these findings indicate that the truncated MLMC estimator performs best when computational resources are limited. 0.675
0
SAA1 SAA2 Untruncated MLMC Truncated MLMC
0.625
−1
0.600
log10(MSE)
Running Average (F̂(x))
0.650
0.575
−2
0.550
Untruncated MLMC Truncated MLMC True value
0.525 0.500 0.475 500
20,000
40,000
60,000
Number of Trees (n1)
80,000
−3
100,000
(a) MLMC estimators with 95% confidence intervals vs n1 .
−4
1
2
3
4
5
log10(Total Sampling Cost)
6
7
(b) Mean squared error (MSE) vs total sampling cost.
Figure 2: Left panel: dependence of the untruncated and truncated MLMC estimators on n1 . Right panel: trade-off between the mean squared error and the total sampling cost of all MLMC and SAA estimators. Each marker corresponds to an average of 10 independent estimators of the same type for a fixed n1 .
Pricing Bermudan Options Consider a Bermudan basket option that grants its holder the right to sell an equally weighted basket of m = 5 assets at a fixed strike price K = 100 on one of T = 4 predetermined exercise dates t ∈ [T ]. We use ξt ∈ Rm to denote the random vector of asset prices at time t and assume that the risk-free interest to γ = 0.05. Exercising the option at time t ∈ [T ] yields an immediate payoff rate amounts 1 ⊤ gt (ξt ) = max 0, K − m em ξt , where em ∈ Rm denotes the vector of ones. Under the risk-neutral measure, the asset prices follow independent geometric Brownian motions with drift γ and volatility σ = 0.2. Hence, the asset prices at the exercise dates satisfy ξt+1 = diag(εt ) ξt , where the disturbances εt ∈ Rm , t ∈ [T − 1], are serially independent and identically distributed with ln(εt ) ∼ N ((γ − σ 2 /2) em , σ 2 Im ). We also assume that ξ1 = 100 em is deterministic. The price of the option is given by the risk-neutral expected net present value of the payoffs under the optimal exercise strategy, which can be found by solving an optimal stopping problem as described in Section 2. If xt is the expected value of the option at time t + 1 conditional on time-t information, then we have xt−1 = Et−1 [max{gt (ξt ), e−γ xt }] for every t ∈ [T ]. Thus, computing today’s price x0 amounts to solving the instance of (MCCE) with ft (ξt , xt ) = max{gt (ξt ), e−γ xt } for all t ∈ [T − 1] and fT (ξT , x) = gT (ξT ). For a more detailed discussion of Bermudan options see [43]. The option pricing problem under consideration does not admit a closed-form solution and must therefore be tackled using sampling-based methods. Approximate prices obtained with classical approaches are reported in [43, Table 4]; for instance, least squares Monte Carlo and the stochastic grid method yield estimates of 2.163(0.001) and 2.141(0.008), respectively, where numbers in parentheses denote standard errors and do not account for bias. In contrast, a policy improvement method [7] produces a 95% confidence in29
terval of [2.154, 2.164] for the true option value. All of these classical methods are specifically designed for option pricing and do not readily extend to general MCCE problems. Furthermore, as noted in [72], their computational cost typically scales superlinearly with the dimension m, which limits their applicability to high-dimensional basket options, where m can be as large as 5,000 [6]. In [72] an untruncated MLMC estimator is shown to produce an unbiased price estimate of 2.161 for the Bermudan option under consideration. Further numerical experiments suggest that the computational cost of this estimator scales sublinearly with m. However, the analysis of this estimator heavily exploits the special structure of optimal stopping problems. As for any MLMC estimator, its efficiency depends critically on the rate parameters rt governing the geometric distributions of the log-branching factors. If rt ≤ 0.5, then the estimator’s expected cost diverges. In addition, there exists a positive tolerance δ ≤ 10−4 such that the estimator’s variance cannot be guaranteed to be finite for any rt ≥ 0.5 + δ [72, Theorem 2]. Within the theoretically admissible narrow range of rt close to 0.5, the untruncated MLMC estimator typically incurs high computational costs. For instance, setting rt = 0.5001 for t ∈ [T − 1] (δ = 10−4 ) results in an expected cost of 2500.5T −1 per tree. We now compute the option price using three truncated MLMC estimators with Mt = 9, 10, 11, as proposed in this paper, and two untruncated MLMC estimators with Mt = ∞ proposed in [72]. The selected truncation points for the truncated estimators empirically provide the best balance between bias and computational cost. For the truncated MLMC estimators, we adopt a heuristic proposed in [72] to choose identical rate parameters rt = r for all geometric distributions governing the log-branching factors. As discussed after Lemma 3.22, the rate r should approximately minimize the expected sampling cost while 2 /2. This is achieved by minimizing the workensuring that the estimator’s variance does not exceed QT −1ϵP Mt ℓ 2 normalized second moment [11], defined as µ1 (x) t=1 ℓ=0 qt (ℓ)2 (see also the proof of Lemma C.7). We approximate this quantity for each r ∈ {0.51, 0.52, . . . , 0.70} by estimating the expectation in the b 1 (x), each computed as definition (19) of µ21 (x) via a sample average over 106 independent copies of H in Definition 3.17 but from a simplified and computationally cheap option pricing problem in which the 1 ⊤ average asset price m e ξt is replaced by a standard normal distribution. Next, we fit a piecewise linear convex function of r to the resulting approximate work-normalized second moments via a least squares procedure and select the rate parameter for all geometric distributions as a minimizer of this function. This yields r = 0.59 for the truncation points Mt = 9 and Mt = 11 and r = 0.58 for Mt = 10. These rates are consistent with the theory of Section 3.2, that is, they yield estimators with finite biases, variances and expected costs. As in [72], we use the same heuristic approach to find a near-optimal rate of r = 0.60 for the first untruncated MLMC estimator. However, this rate is inconsistent with the theory of [72], that is, for this choice of r the variance of the estimator is not guaranteed to be finite. The rate of the second untruncated MLMC estimator is set to r = 0.5001, which is the largest rate that guarantees the estimator to have a finite variance. Table 1 reports the option prices predicted by all truncated and untruncated MLMC estimators when the number of scenario trees is set to n1 = 5 × 106 . Even though they only use half as many scenario trees as the methods in [7, 72], our truncated MLMC estimators provide comparable 95% confidence intervals. Note that the untruncated MLMC estimator with r = 0.6 offers no confidence interval because its variance is not guaranteed to be finite, while the untruncated MLMC estimator with r = 0.5001 suffers from an excessive sampling cost per tree. Contextual Bandits In the last experiment we solve an instance of the off-policy learning problem (12) borrowed from [59], which determines one of two treatments to be administered to hospital patients suffering from acute ischaemic stroke. Each patient is characterized by a context vector c ∈ C ⊆ Z6+ . The first context feature c1 ∈ {0, 1, 2} represents the patient’s level of consciousness (0: fully conscious, 1: drowsy, 2: unconscious), while c2 ∈ {0, . . . , 4}, c3 , c4 ∈ {0, 1}, c5 ∈ {0, . . . , 5} and c6 ∈ {0, . . . , 3} capture other categorical features such as the patient’s age etc. The overall context space C is thus finite, and its cardinality equals 1,440. As in [59], we assume that there are two actions a ∈ A = {1, 2} (1: prescribing both aspirin 30
Mt 9 10 11 ∞ ∞
rt 0.5900 0.5800 0.5900 0.5001 0.6000
Cost per tree 22.6084 29.5795 26.3283 1.5634 × 1010 27.0000
Estimate (se) 2.1684(0.0076) 2.1562(0.0072) 2.1641(0.0080) 2.1547(0.0073) 2.1589(−)
Confidence interval [2.1536, 2.1832] [2.1421, 2.1703] [2.1485, 2.1797] [2.1403, 2.1690] −
Table 1: MLMC estimators for pricing Bermudan basket options: expected costs per tree, estimated option prices with standard errors (se) and 95% confidence intervals for different truncation points and rates. and heparin, 2: not administering any treatment), and we parametrize the context-dependent policy πθ with θ ∈ [0, 1]2 , where the agent selects action 1 with probability θ1 if c1 ̸= 0 and with probability θ2 if c1 = 0. We will solve the resulting instance of (MCCO) with a variant of the Adam optimizer [45] equipped with the MLMC gradient estimator of Definition 4.2. This allows us to assess our MLMC gradient estimator in an interesting scenario where we have no theoretical results guaranteeing its superiority over basic SAA estimators. Indeed, as explained in Section 2, the integrands f1 , f2 and f3 involve logarithms as well as quadratic and exponential terms. Thus, they fail to satisfy the Lipschitz continuity and smoothness conditions detailed in Assumptions 3.2 and 3.3, which are needed for all theoretical results in Section 4. We prescribe the joint distribution of c, u and y under P synthetically, which allows us to solve problem (12) exactly. The resulting ground-truth solution serves as a benchmark against which any data-driven policies can be compared. First, we assume that the context c as well as the auxiliary random vector u follow independent uniform distributions on C. Next, the mean of the cost vector y conditional on c is set to (3.0 + 5.0c5 + p(c), 5.5 + 1.0c5 + p(c)) if c1 = 0, E[y|c] = (1.7 + 3.5c5 + p(c), 3.0 + 1.0c5 + p(c)) if c1 ̸= 0, where p(c) = 2.4+1.92(c5 /5−2.5)2 if c2 = 4, c3 = 1, c4 = 1 and c6 = 3; p(c) = 0 otherwise. Conditional on c, the actual cost vector y follows a bivariate log-normal distribution, that is, log(y) follows the normal distribution with mean log(E[y|c])− 12 diag(Σ) and covariance matrix Σ, where the diagonal and off-diagonal elements of Σ are set to 5.0 and 2.5, respectively. The heterogeneity of E[y|c] ensures that the exact optimal solution of problem (12) displays a sufficient level of complexity, whereas the positive correlation between the costs of the two actions captures the effect of an unobserved confounder that simultaneously affects both treatment outcomes. Finally, we set the smoothness temperature in (12) to µ = 2, and we set the ambiguity radii to rc = 0.4 and ry = 0.15. This choice roughly accounts for a covariate shift whereby the subpopulation of patients with c1 = 0 is over-sampled by 50%. By leveraging known properties of log-normal distributions, the innermost integral in (12) can be evaluated analytically. As C is finite, the objective of (12) thus reduces to a standard log-sum-exp function. Hence, the exact global minimizer (λ∗ , θ1∗ , θ2∗ ) ≈ (11.829, 0.589, 0.713) of problem (12) can be computed efficiently to any desired accuracy with any off-the-shelf convex optimization solver. We now address problem (12) with oracle-based methods that have only sample access to P. Specifically, we solve (12) with a variant of the Adam optimizer using learning rates of 0.75 for λ and 0.025 for θ1 and θ2 , alongside default decay rate parameters. To prevent Adam’s momentum from accumulating against a hard boundary and stalling convergence, we enforce the λ ≥ 0 constraint via a softplus reparameterization rather than a standard projection. Concretely, we set λ = log(1 + exp(λ′ )) for some λ′ ∈ R and run Adam entirely on λ′ . However, we simply project θ1 and θ2 to [0, 1] in each iteration of the algorithm. We also add the L2 regularization term 0.005·∥θ∥22 to the objective of (12) to enhance stability. As stochastic gradient estimators we use the SAA-type estimator proposed in [59] (adapted to account for stochastic costs) as well as the MLMC gradient estimator of Definition 4.2. As the Lipschitz continuity and smoothness conditions re31
quired to control the bias and variance of the MLMC estimator fail to hold in this example, we clip the norms of the truncated MLMC gradient estimators with respect to λ, θ1 and θ2 at 100, 50 and 50, respectively. Note that problem (12) as described in [59] constitutes an instance of (MCCO) with T = 3. However, the numerical experiments in [59] focus on a simplified model with deterministic costs and T = 2—presumably due to the hardness of the original three-stage problem. In contrast, our experiments faithfully capture the three-stage nature of the more challenging original problem with stochastic costs. Figure 3 illustrates the convergence of λ, θ1 and θ2 as a function of the cumulative number of scenarios over 2,000 Adam iterations and for different choices of the estimators’ hyperparameters. Solid lines and shaded regions represent means as well as corresponding 95% confidence intervals obtained from 20 independent simulation runs, whereas dotted lines represent ground-truth minimizers. We observe that Adam converges significantly faster when equipped with MLMC estimators instead of SAA estimators. Recall from Lemma 4.6 that the truncated MLMC estimator is biased, and note that gradient clipping introduces an additional bias. Nevertheless, MLMC estimators lead to less biased optimizers than SAA estimators. Even though Assumptions 3.2 and 3.3 fail to hold, MLMC estimators thus enable faster and tighter convergence to the true optimal solution than SAA estimators equipped with a comparable computational budget. SAA Setups (n1 , n2 , n3 ) (100,200,200)
60
(100,300,300)
MLMC Setups (n1 , M, r)
(100,500,500)
(100,1000,1000)
(1000,(9,5),(0.5,0.5))
0.8
True Optimum λ (11.829)
(1000,(10,6),(0.5,0.5))
True Optimum θ1 (0.589)
True Optimum θ2 (0.713)
0.7
50
0.6 40 30
0.4
θ2
λ
θ1
0.5
0.3 20
0.2 0.1
10 5 105
106
107
108
# Scenarios
109
1010
1011
0.0 105
106
107
108
109
# Scenarios
1010
1011
105
106
107
108
109
# Scenarios
1010
1011
Figure 3: Convergence of the decision variables λ (left), θ1 (middle), and θ2 (right) generated by Adam. Results show mean trajectories with 95% confidence intervals based on 20 independent simulation runs.
Acknowledgments This work was supported as a part of NCCR Automation, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 51NF40 225155). We thank Guanyang Wang for providing the code for the experiments in [61], as well as Yi Shen and Pan Xu for the code corresponding to [59]. We also thank Guanyang Wang, Lie He and Niao He for helpful discussions.
References [1] Charalambos D Aliprantis and Kim C Border. Infinite Dimensional Analysis: A Hitchhiker’s Guide. Springer, 2006. [2] Yossi Arjevani, Yair Carmon, John C Duchi, Dylan J Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, 199(1–2):165– 214, 2023.
32
[3] Waı̈ss Azizian, Franck Iutzeler, and Jérôme Malick. Regularization for Wasserstein distributionally robust optimization. ESAIM: Control, Optimisation and Calculus of Variations, 29:1–33, 2023. [4] Krishnakumar Balasubramanian, Saeed Ghadimi, and Anthony Nguyen. Stochastic multilevel composition optimization algorithms with level-independent convergence rates. SIAM Journal on Optimization, 32(2):519–544, 2022. [5] Joakim Beck, Ben Mansour Dia, Luis Espath, and Raúl Tempone. Multilevel double loop Monte Carlo and stochastic collocation methods with importance sampling for Bayesian optimal experimental design. International Journal for Numerical Methods in Engineering, 121(15):3482–3503, 2020. [6] Sebastian Becker, Patrick Cheridito, Arnulf Jentzen, and Timo Welti. Solving high-dimensional optimal stopping problems using deep learning. European Journal of Applied Mathematics, 32(3):470– 514, 2021. [7] Christian Bender, Anastasia Kolodko, and John Schoenmakers. Policy iteration for American options: Overview. Monte Carlo Methods and Applications, 12(5):347–362, 2006. [8] Andrew Bennett, Nathan Kallus, and Tobias Schnabel. Deep generalized method of moments for instrumental variable analysis. In Advances in Neural Information Processing Systems, pages 3564– 3574, 2019. [9] Dimitri P Bertsekas. Dynamic Programming and Optimal Control, volume 1. Athena Scientific, 3rd edition, 1995. [10] Jose Blanchet, Donald Goldfarb, Garud Iyengar, Fengpei Li, and Chaoxu Zhou. Unbiased simulation for optimizing stochastic function compositions. arXiv:1711.07564, 2017. [11] Jose H Blanchet and Peter W Glynn. Unbiased Monte Carlo for optimization and functions of expectations via multi-level randomization. In Winter Simulation Conference, pages 3656–3667, 2015. [12] Mark Broadie, Yiping Du, and Ciamac C Moallemi. Efficient risk estimation via nested sequential simulation. Management Science, 57(6):1172–1194, 2011. [13] Karolina Bujok, Ben M Hambly, and Christoph Reisinger. Multilevel simulation of functionals of Bernoulli random variables with application to basket credit derivatives. Methodology and Computing in Applied Probability, 17:579–604, 2015. [14] Tianyi Chen, Yuejiao Sun, and Wotao Yin. Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization. IEEE Transactions on Signal Processing, 69:4937–4948, 2021. [15] Xuxing Chen, Abhishek Roy, Yifan Hu, and Krishnakumar Balasubramanian. Stochastic optimization algorithms for instrumental variable regression with streaming data. In Advances in Neural Information Processing Systems, pages 26510–26542, 2024. [16] Dragos Florin Ciocan and Velibor V Mišić. Interpretable optimal stopping. Management Science, 68 (3):1616–1638, 2022. [17] Weilin Cong, Rana Forsati, Mahmut Kandemir, and Mehrdad Mahdavi. Minimal variance sampling with provable guarantees for fast training of graph neural networks. In International Conference on Knowledge Discovery & Data Mining, pages 1393–1403, 2020. 33
[18] Bo Dai, Niao He, Yunpeng Pan, Byron Boots, and Le Song. Learning from conditional distributions via dual embeddings. In Artificial Intelligence and Statistics, pages 1458–1467, 2017. [19] Bo Dai, Albert Shaw, Lihong Li, Lin Xiao, Niao He, Zhen Liu, Jianshu Chen, and Le Song. SBEED: Convergent reinforcement learning with nonlinear function approximation. In International Conference on Machine Learning, pages 1125–1134, 2018. [20] Martin Dyer and Leen Stougie. Computational complexity of stochastic programming problems. Mathematical Programming, 106(3):423–432, 2006. [21] Hongchang Gao. Decentralized multi-level compositional optimization algorithms with levelindependent convergence rate. In International Conference on Artificial Intelligence and Statistics, pages 4402–4410, 2024. [22] Saeed Ghadimi, Guanghui Lan, and Hongchao Zhang. Mini-batch stochastic approximation methods for nonconvex stochastic composite optimization. Mathematical Programming, 155(1):267–305, 2016. [23] Michael B Giles. Multilevel Monte Carlo path simulation. Operations Research, 56(3):607–617, 2008. [24] Michael B Giles. Multilevel Monte Carlo methods. Acta Numerica, 24:259–328, 2015. [25] Michael B Giles. MLMC for nested expectations. In Josef Dick, Frances Y. Kuo, and Henryk Woźniakowski, editors, Contemporary Computational Mathematics: A Celebration of the 80th Birthday of Ian Sloan, pages 425–442. Springer, 2018. [26] Michael B Giles and Abdul-Lateef Haji-Ali. Multilevel nested simulation for efficient risk estimation. SIAM/ASA Journal on Uncertainty Quantification, 7(2):497–525, 2019. [27] Michael B Giles and Lukasz Szpruch. Antithetic multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation. Annals of Applied Probability, 24(4):1585–1620, 2014. [28] Michael B Giles, Abdul-Lateef Haji-Ali, and Jonathan Spence. Efficient risk estimation for the credit valuation adjustment. arXiv:2301.05886, 2023. [29] Takashi Goda and Wataru Kitade. Constructing unbiased gradient estimators with finite variance for conditional stochastic optimization. Mathematics and Computers in Simulation, 204:743–763, 2023. [30] Takashi Goda, Tomohiko Hironaka, and Takeru Iwamoto. Multilevel Monte Carlo estimation of expected information gains. Stochastic Analysis and Applications, 38(4):581–600, 2020. [31] Takashi Goda, Tomohiko Hironaka, Wataru Kitade, and Adam Foster. Unbiased MLMC stochastic gradient-based optimization of Bayesian experimental designs. SIAM Journal on Scientific Computing, 44(1):A286–A311, 2022. [32] Michael B Gordy and Sandeep Juneja. Nested simulation in portfolio risk measurement. Management Science, 56(10):1833–1848, 2010. [33] Abdul-Lateef Haji-Ali and Jonathan Spence. Nested multilevel Monte Carlo with biased and antithetic sampling. arXiv:2308.07835, 2023. [34] Grani A Hanasusanto, Daniel Kuhn, and Wolfram Wiesemann. A comment on “computational complexity of stochastic programming problems”. Mathematical Programming, 159(1–2):557–569, 2016. 34
[35] Lars Peter Hansen and Thomas J Sargent. Robustness. Princeton University Press, 2008. [36] Jason Hartford, Greg Lewis, Kevin Leyton-Brown, and Matt Taddy. Deep IV: A flexible approach for counterfactual prediction. In International Conference on Machine Learning, pages 1414–1423, 2017. [37] Lie He and Shiva Kasiviswanathan. Debiasing conditional stochastic optimization. In Advances in Neural Information Processing Systems, pages 78846–78893, 2023. [38] Roger A Horn and Charles R Johnson. Matrix Analysis. Cambridge University Press, 1985. [39] Yifan Hu, Xin Chen, and Niao He. Sample complexity of sample average approximation for conditional stochastic optimization. SIAM Journal on Optimization, 30(3):2103–2133, 2020. [40] Yifan Hu, Siqi Zhang, Xin Chen, and Niao He. Biased stochastic first-order methods for conditional stochastic optimization and applications in meta learning. In Advances in Neural Information Processing Systems, pages 2759–2770, 2020. [41] Yifan Hu, Xin Chen, and Niao He. On the bias-variance-cost tradeoff of stochastic optimization. In Advances in Neural Information Processing Systems, pages 22119–22131, 2021. [42] Yifan Hu, Jie Wang, Xin Chen, and Niao He. Multi-level Monte-Carlo gradient methods for stochastic optimization with biased oracles. arXiv:2408.11084, 2024. [43] Shashi Jain and Cornelis W. Oosterlee. Pricing high-dimensional Bermudan options using the stochastic grid method. International Journal of Computer Mathematics, 89(9):1186–1211, 2012. [44] Wei Jiang, Bokun Wang, Yibo Wang, Lijun Zhang, and Tianbao Yang. Optimal algorithms for stochastic multi-level compositional optimization. In International Conference on Machine Learning, pages 10195–10216, 2022. [45] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. [46] Daniel Kuhn, Soroosh Shafiee, and Wolfram Wiesemann. Distributionally robust optimization. Acta Numerica, 34:579–804, 2025. [47] Tze Leung Lai. Optimal stopping and sequential tests which minimize the maximum expected sample size. Annals of Statistics, pages 659–673, 1973. [48] Yifan Lin, Yuxuan Ren, and Enlu Zhou. Bayesian risk Markov decision processes. In Advances in Neural Information Processing Systems, pages 17430–17442, 2022. [49] Krikamol Muandet, Arash Mehrjou, Si Kai Lee, and Anant Raj. Dual instrumental variable regression. In Advances in Neural Information Processing Systems, pages 2710–2721, 2020. [50] Emin Ozyoruk, Nesim Kohen Erkip, and Çağın Ararat. End-of-life inventory management problem: Results and insights. International Journal of Production Economics, 243:108313, 2022. [51] Tom Rainforth, Rob Cornish, Hongseok Yang, Andrew Warrington, and Frank Wood. On nesting Monte Carlo estimators. In International Conference on Machine Learning, pages 4267–4276, 2018. [52] Marcus de Mendes C. R. Reaiche. A note on sample complexity of multistage stochastic programs. Operations Research Letters, 44(4):430–435, 2016. 35
[53] Chang-han Rhee and Peter W Glynn. Unbiased estimation with square root convergence for SDE models. Operations Research, 63(5):1026–1043, 2015. [54] Andrzej Ruszczynski. A stochastic subgradient method for nonsmooth nonconvex multilevel composition optimization. SIAM Journal on Control and Optimization, 59(3):2301–2320, 2021. [55] Andrzej Ruszczyński and Alexander Shapiro. Conditional risk mappings. Mathematics of Operations Research, 31(3):544–561, 2006. [56] Alexander Shapiro. On complexity of multistage stochastic programs. Operations Research Letters, 34(1):1–8, 2006. [57] Alexander Shapiro and Arkadi Nemirovski. On complexity of stochastic programming problems. In Vaithilingam Jeyakumar and Alexander Rubinov, editors, Continuous Optimization, pages 111–146. Springer, Boston, MA, 2005. [58] Alexander Shapiro, Enlu Zhou, and Yifan Lin. Bayesian distributionally robust optimization. SIAM Journal on Optimization, 33(2):1279–1304, 2023. [59] Yi Shen, Pan Xu, and Michael M Zavlanos. Wasserstein distributionally robust policy evaluation and learning for contextual bandits. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. Featured Certification. [60] Rahul Singh, Maneesh Sahani, and Arthur Gretton. Kernel instrumental variable regression. In Advances in Neural Information Processing Systems, pages 4593–4605, 2019. [61] Yasa Syed and Guanyang Wang. Optimal randomized multilevel Monte Carlo for repeatedly nested expectations. In International Conference on Machine Learning, pages 33343–33364, 2023. [62] Emanuel Todorov and Michael I. Jordan. Optimal feedback control as a theory of motor coordination. Nature Neuroscience, 5(11):1226–1235, 2002. [63] Nicky D Van Foreest and Onur A Kilic. An intuitive approach to inventory control with optimal stopping. European Journal of Operational Research, 311(3):921–924, 2023. [64] Martin J. Wainwright. High-Dimensional Statistics: A Non-Asymptotic Viewpoint. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2019. [65] Jie Wang, Rui Gao, and Yao Xie. Sinkhorn distributionally robust optimization. Operations Research, 2025. Forthcoming. [66] Tianze Wang and Guanyang Wang. Unbiased Multilevel Monte Carlo methods for intractable distributions: MLMC meets MCMC. Journal of Machine Learning Research, 24(249):1–40, 2023. [67] Yuhao Wang and Enlu Zhou. Bayesian risk-averse Q-learning with streaming observations. In Advances in Neural Information Processing Systems, pages 75967–75992, 2024. [68] Peter Whittle. Risk-Sensitive Optimal Control. Wiley, 1990. [69] Tesi Xiao, Krishnakumar Balasubramanian, and Saeed Ghadimi. A projection-free algorithm for constrained stochastic multi-level composition optimization. In Advances in Neural Information Processing Systems, pages 19984–19996, 2022.
36
[70] Shuoguang Yang, Mengdi Wang, and Ethan X Fang. Multilevel stochastic gradient methods for nested composition optimization. SIAM Journal on Optimization, 29(1):616–659, 2019. [71] Junyu Zhang and Lin Xiao. Multilevel composite stochastic optimization via nested variance reduction. SIAM Journal on Optimization, 31(2):1131–1157, 2021. [72] Zhengqing Zhou, Guanyang Wang, Jose H Blanchet, and Peter W Glynn. Unbiased optimal stopping via the MUSE. Stochastic Processes and their Applications, 166:104088, 2023.
Appendix A
Auxiliary Results
The following lemma establishes a uniform deviation bound based on covering numbers. It is a standard result in stochastic programming, and we include a concise proof to keep the presentation self-contained. Lemma A.1 ([56, Proof of Theorem 1]). Consider the problem of minimizing an L-Lipschitz continuous function F : Rd → R over a compact set X ⊆ Rd with diameter DX . Let Fb : Rd → R be an L-Lipschitz continuous random function that approximates F . Then, for every v > 0, there is a v-net Xv ⊆ X such that ϵ |Xv | ≤ ⌈2DX /v + 1⌉d and minx′ ∈Xv ∥x − x′ ∥2 ≤ v for every x ∈ X . If ϵ > 0 and v = 4L , then we have X ϵ P sup |Fb(x) − F (x)| > ϵ ≤ P |Fb(x′ ) − F (x′ )| > . 2 x∈X ′ x ∈Xv
Proof. The set X is contained in a ball of radius DX and can thus be covered by ⌈2DX /v + 1⌉d smaller balls of diameter v [64, Lemma 5.7]. By construction, the centroids of these smaller balls, projected onto X , constitute a v-net Xv . That is, for every x ∈ X , there is an x′ ∈ Xv in its closed v-neighborhood. Assume ϵ now that v = 4L for some ϵ > 0. Then, for every x ∈ X and x′ ∈ argminx′ ∈Xv ∥x − x′ ∥2 , we have ϵ |Fb(x) − F (x)| ≤ |Fb(x) − Fb(x′ )| + |Fb(x′ ) − F (x′ )| + |F (x′ ) − F (x)| ≤ |Fb(x′ ) − F (x′ )| + , 2 where the second inequality follows from the L-Lipschitz continuity of F and Fb and from the observation that ∥x′ − x∥2 ≤ v. As the above bound holds for every x ∈ X , we readily obtain sup |Fb(x) − F (x)| ≤
x∈X
ϵ + max |Fb(x′ ) − F (x′ )|. 2 x′ ∈Xv
This further implies via the union bound that X ϵ ′ ′ b(x ) − F (x )| > b(x′ ) − F (x′ )| > ϵ . P sup |Fb(x) − F (x)| > ϵ ≤ P max | F ≤ P | F x′ ∈Xv 2 2 x∈X ′ x ∈Xv
Hence, the claim follows. The following matrix norm bounds will be used repeatedly in proving the main results of this paper. 1
Lemma A.2 ([38, § 5.4.P3]). If A ∈ Rm×n and 1 ≤ p < q ≤ ∞, then ∥A∥q ≤ ∥A∥p ≤ (mn) p Lemma A.3 ([38, § 5.6.0.2]). If A ∈ Rm×n and B ∈ Rn×k , then ∥AB∥2 ≤ ∥A∥2 ∥B∥2 . 37
− 1q
∥A∥q .
Finally, we state a well-known corollary of the Marcinkiewicz-Zygmund inequality as an auxiliary lemma, which is useful for deriving variance bounds of the estimators proposed in the main text. Lemma A.4. Assume that Y and Z are random matrices of arbitrary dimensions with finite p-th moments for some p ≥ 2, and let Zi , i ∈ [n], be independent samples from the conditional distribution PZ|Y . Then there exists a constant Bp depending only on p such that EZ|Y
n p Bp 1X Zi − EZ|Y [Z] ≤ p/2 EZ|Y ∥Z∥pp . n p n i=1
The proof parallels that of [72, Corollary 2], generalized to matrices. Details are omitted for brevity.
B
Proofs of Section 2
Proof of Lemma 2.1. We prove the claim by backward induction on t. At t = T we have QT (s, a) = s⊤ PT s by assumption, which is of the form (6a) if the coefficients satisfy (6c). This establishes the base step. As for the induction step, assume that (6a) holds for Qt+1 . Thus, we find ⊤
⊤ s′ + dt+1 , min Qt+1 (s′ , a′ ) = s′ Pt+1 s′ + 2 gt+1
a′ ∈Rn
with the Schur-complement expressions for Pt+1 , gt+1 and dt+1 given in (6d). Substituting the above formula with s′ = As + Ba + ξt+1 into the Bellman recursion (4) and expanding the quadratic and linear terms shows that Qt (s, a) is again a quadratic function in (s, a). Collecting coefficients of s⊤ s, s⊤ a, a⊤ a, s, a and the constant term then yields exactly the expressions (6b). Thus the quadratic structure of the Q-function is preserved backward in time, completing the induction. The proof of Lemma 2.1 readily reveals that the optimal control at time t is given by a⋆t = − Qaa t
−1
⊤ a (Qsa t ) s + bt
and thus constitutes an affine function of the state s with Ft -measurable random coefficients.
C
Proofs of Section 3
To streamline the proofs in Section 3, we first introduce some key notation and conventions that will be useful for our derivations. Throughout this discussion we use Lt as shorthand for the space L1 (Ω, F, P; Rdt ) of all P-integrable dt -dimensional random vectors on (Ω, F, P). Definition C.1. For any t ∈ [T ] and s ∈ [t + 1], we define the conditional risk mapping Fst : Lt → Ls−1 via ( Es−1 fs ξs , . . . , Et−1 ft (ξt , zt ) . . . if s < t + 1, t Fs (zt ) := zt if s = t + 1. Note that Fst maps dt -dimensional random vectors to ds−1 -dimensional random vectors. By construction, we have F (x) = F1T (x) for every (deterministic) vector x ∈ Rd . If Assumptions 3.1 and 3.2 hold, then Fst (zt ) is well-defined and constitutes indeed an element of Ls−1 for every zt ∈ Lt . Lemma C.2. If Assumptions 3.1 and 3.2 hold, then Fst (zt ) ∈ Ls−1 for all zt ∈ Lt , t ∈ [T ] and s ∈ [t + 1]. 38
Proof. We fix t ∈ [T ] and zt ∈ Lt and prove the claim by induction on s. If s = t + 1, then we have Fst (zt ) = zt ∈ Lt , and thus the base case is trivially true. As for the induction step corresponding to t (z ) ∈ L . Then, f (ξ , F t (z )) is F-measurable thanks to [1, Lemma 4.51] any s ∈ [t], assume that Fs+1 t s s s s+1 t and because measurability is preserved under Borel-measurable transformations. Also, we have t E Fst (zt ) 2 = E Es−1 [fs (ξs , Fs+1 (zt ))] 2 t ≤ E fs (ξs , Fs+1 (zt )) 2 t (zt )) − fs (ξs , 0) 2 + fs (ξs , 0) 2 ≤ E fs (ξs , Fs+1 t (zt ) 2 + E fs (ξs , 0) 2 < ∞, ≤ Ls E Fs+1 where the first inequality follows from Jensen’s inequality and the law of total expectation, while the second and the third inequalities follow from the triangle inequality and Assumption 3.2, respectively. The resulting t (z ) ∈ L , and by Assumption 3.1. We bound is finite by the induction hypothesis, which assumes that Fs+1 t s t may thus conclude that Fs (zt ) ∈ Ls−1 . This observation completes the induction step. Lemma C.2 ensures that the function F (x) is well-defined. If Assumptions 3.1 and 3.2 hold, then all conditional risk mappings of Definition C.1 are Lipschitz continuous. If additionally Assumption 3.3 holds, then the conditional risk mappings are also Lipschitz smooth. Lemma C.3. Fix any t ∈ [T ] and s ∈ [t + 1]. If Assumptions 3.1 and 3.2 hold, then Fst is an L[s:t] -Lipschitz Q continuous function on the space Lt with L[s:t] := ti=s Li if s < t + 1; := 1 if s = t + 1. This means that Fst (zt ) − Fst (zt′ ) 2 ≤ L[s:t] Es−1 ∥zt − zt′ ∥2
∀zt , zt′ ∈ Lt .
If additionally Assumption 3.3 holds, then Fst is an S[s:t] -smooth function on the space Lt with S[s:t] := Pt 2 i=s Si L[s:i−1] L[i+1:t] if s < t + 1; := 0 if s = t + 1. This means that ∇Fst (zt ) − ∇Fst (zt′ ) 2 ≤ S[s:t] Es−1 ∥zt − zt′ ∥2
∀zt , zt′ ∈ Lt .
Proof. We fix t ∈ [T ] and use backward induction on s to prove that Fst is L[s:t] -Lipschitz continuous. In the base case corresponding to s = t + 1, Fst collapses to the identity mapping on Lt , and thus the claim is t trivially true. As for the induction step corresponding to any s ∈ [t], assume that Fs+1 is L[s+1:t] -Lipschitz t t continuous. Since Fs (zt ) = Es−1 [fs (ξs , Fs+1 (zt ))], we thus have t t Fst (zt ) − Fst (zt′ ) 2 ≤Es−1 fs ξs , Fs+1 (zt ) − fs ξs , Fs+1 (zt′ ) 2 t t ≤Ls Es−1 Fs+1 (zt ) − Fs+1 (zt′ ) 2 ≤Ls L[s+1:t] Es−1 ∥zt − zt′ ∥2 = L[s:t] Es−1 ∥zt − zt′ ∥2 for any zt , zt′ ∈ Lt . Here, the first two inequalities follow from Jensen’s inequality and from Assumption 3.2, respectively, while the third inequality exploits the induction hypothesis and the law of total expectation. This completes the induction step, and thus the claim follows. Suppose now that Assumption 3.3 holds, too. One readily verifies that the smoothness constants obey the recursion S[s:t] = Ls S[s+1:t] + Ss L2[s+1:t] for all s ∈ [t] and t ∈ [T ]. In the remainder we fix any t ∈ [T ] and use backward induction on s to prove that Fst is S[s:t] -smooth. In the base case corresponding to s = t + 1, the claim is trivially true because Fst is the identity mapping on Lt . As for the induction step corresponding t to any s ∈ [t], assume that Fs+1 is S[s+1:t] -smooth. By the definition of Fst (zt ), we have t t ∇Fst (zt ) = ∇Es−1 fs ξs , Fs+1 (zt ) = Es−1 ∇zt fs ξs , Fs+1 (zt ) . 39
Note that the gradient may be interchanged with the conditional expectation by the dominated convergence theorem, which applies thanks to Assumption 3.2. Jensen’s inequality and the chain rule then imply that t t (zt′ ) 2 (zt ) − ∇zt fs ξs , Fs+1 ∇Fst (zt ) − ∇Fst (zt′ ) 2 ≤ Es−1 ∇zt fs ξs , Fs+1 t t t t (zt′ ) 2 . (zt ) − ∇Fs+1 (zt′ )∇xs fs ξs , Fs+1 = Es−1 ∇Fs+1 (zt )∇xs fs ξs , Fs+1 Next, we may use the triangle inequality to conclude that t t t (zt ) 2 ∇Fst (zt ) − ∇Fst (zt′ ) 2 ≤ Es−1 ∇Fs+1 (zt ) − ∇Fs+1 (zt′ ) ∇xs fs ξs , Fs+1 t t t (zt′ ) (zt ) − ∇xs fs ξs , Fs+1 + Es−1 ∇Fs+1 (zt′ ) ∇xs fs ξs , Fs+1 t t ≤ Ls Es−1 ∇Fs+1 (zt ) − ∇Fs+1 (zt′ ) 2 t t t + Ss Es−1 ∇Fs+1 (zt′ ) 2 Fs+1 (zt ) − Fs+1 (zt′ ) 2 ≤ Ls S[s+1:t] + Ss L2[s+1:t] zt − zt′ 2 = S[s:t] zt − zt′ 2 ,
2
where the second inequality follows from Assumptions 3.2 and 3.3 and from the Cauchy-Schwarz inequality. t The third inequality exploits the induction hypothesis, whereby Fs+1 is S[s+1:t] -smooth, as well as the first t part of the proof, whereby Fs+1 is L[s+1:t] -Lipschitz continuous. Finally, the equality follows from the recursive formula for S[s:t] established earlier. Hence, Fst (zt ) is indeed S[s:t] -smooth in zt . In analogy to Definition C.1, we can now introduce empirical conditional risk mappings on Lt .
Definition C.4. For any t ∈ [T ] and s ∈ [t + 1], we define the conditional risk mapping Fbst : Lt → Ls−1 via ( b s−1 fs ξs , . . . , E b t−1 ft (ξt , zt ) . . . E if s < t + 1, t b Fs (zt ) := zt if s = t + 1. b t−1 [·], t ∈ [T ], denote the empirical conditional expectation operators defined in (15). Here, E
Note that Fb(x) = Fb1T (x) for every (deterministic) vector x ∈ Rd . If Assumptions 3.1 and 3.2 hold, one can proceed as in Lemmas C.2 and C.3 to show that the conditional risk mapping Fbst is well-defined and L[s:t] -Lipschitz continuous for any t ∈ [T ] and s ∈ [t+1]. In fact, the proofs simplify because P-integrability and Lipschitz continuity are trivially preserved under finite sums. Details are omitted for brevity. We are now ready to prove the formal results of Section 3.1. Proof of Lemma 3.6. Suppose first that only Assumptions 3.1, 3.2 and 3.5 hold. By the construction of the conditional risk mappings in Definitions C.1 and C.4, we have F (x) = F1T (x) = F1T (FbTT+1 (x))
and
E[Fb(x)] = E[Fb1T (x)] = F11 (Fb2T (x))
for all x ∈ Rd . We can thus represent the bias of the SAA estimator Fb(x) as a telescoping sum and then use the triangle inequality to conclude that E[Fb(x)] − F (x) ≤
T −1 X t=1
T T F1t Fbt+1 (x) − F1t+1 Fbt+2 (x) .
The t-th term in the resulting upper bound satisfies T T T T F1t Fbt+1 (x) − F1t+1 Fbt+2 (x) = F1t Fbt+1 (x) − F1t Et [Fbt+1 (x)] T T ≤L[t] E Fbt+1 (x) − Et [Fbt+1 (x)] 2 T 2 1 T ≤L[t] E Fbt+1 (x) − Et [Fbt+1 (x)] 2 2 , 40
(27)
where the two inequalities follow from the Lipschitz continuity of F1t established in Lemma C.3 and from Jensen’s inequality, respectively. By the definition of FbtT , we then have T 2 1 T L[t] E Fbt+1 (x) − Et [Fbt+1 (x)] 2 2 h 2 i 12 b t ft+1 ξt+1 , FbT (x) − Et ft+1 ξt+1 , FbT (x) =L[t] E E t+2 t+2 2 h 2 i 21 L[t] σt+1 L[t] T T (x) − Et ft+1 ξt+1 , Fbt+2 (x) E ft+1 ξt+1 , Fbt+2 , ≤ √ =√ nt+1 nt+1 2 where the second equality holds because, conditional on ξ[t] , the variance of an average of nt+1 i.i.d. random variables coincides with the variance of a single random variable divided by nt+1 . The inequality in the last line follows from Assumption 3.5. Substituting the resulting estimate into (27) yields E[Fb(x)] − F (x)| ≤
T −1 X t=1
L[t] σt+1 . √ nt+1
Suppose now that Assumption 3.3 holds, too. Hence, the t-th term in (27) admits the tighter estimate T T (x) F1t Fbt+1 (x) − F1t+1 Fbt+2 T T = F1t−1 Et−1 ft ξt , Fbt+1 (x) − F1t−1 Et−1 ft ξt , Et [Fbt+1 (x)] h h i i T T ≤L[t−1] E Et−1 ft ξt , Fbt+1 (x) − ft ξt , Et [Fbt+1 (x)]
(28)
2
where the inequality follows from the Lipschitz continuity of F1t−1 established in Lemma C.3. h i ⊤ T T T Et−1 ∇xt ft ξt , Et [Fbt+1 Fbt+1 (x) − Et [Fbt+1 (x)] = 0 (x)] T (x)]) is F -measurable. The conditional expectation in the last line of (28) thus satisfies as ∇xt ft (ξt , Et [Fbt+1 t h i T T Et−1 ft ξt , Fbt+1 (x) − ft ξt , Et [Fbt+1 (x)] 2 h i T T T T T b b = Et−1 ft ξt , Ft+1 (x) − ft ξt , Et [Ft+1 (x)] + ∇xt ft (ξt , Et [Fbt+1 (x)])⊤ Fbt+1 (x) − Et [Fbt+1 (x)] 2 h 2i St T T ≤ Et−1 Fbt+1 (x) − Et [Fbt+1 (x)] 2 2 h h i h 2 i 2 i St σt+1 St b t ft+1 ξt+1 , FbT (x) − Et ft+1 ξt+1 , FbT (x) = Et−1 E ≤ . t+2 t+2 2 2nt+1 2
Here, the first inequality follows from the bound in (13), which applies thanks to Assumption 3.3. The second inequality holds again because, conditional on ξ[t] , the variance of an average of nt+1 i.i.d. random variables coincides with the variance of a single random variable divided by nt+1 . Substituting the above estimate into (28) and combining the resulting bound with (27) finally yields E[Fb(x)] − F (x)| ≤ This observation completes the proof. 41
T −1 X t=1
2 L[t−1] St σt+1 . 2nt+1
Proof of Lemma 3.7. Since F (x) is deterministic, the variance of the SAA estimator Fb(x) satisfies h 2 i V Fb(x) = V Fb(x) − F (x) ≤ E Fb(x) − F (x) . To bound the variance of Fb(x), it thus suffices to bound the second moment of the estimation error Fb(x) − F (x). By the construction of the conditional risk mappings in Definitions C.1 and C.4, we have F (x) = F1T (x) = Fb10 (F1T (x))
and Fb(x) = Fb1T (x) = Fb1T (FTT+1 (x))
for all x ∈ Rd . We can thus represent the estimation error as the telescoping sum Fb(x) − F (x) =
T −1 X t=0
T T Fb1t (Ft+1 (x)) − Fb1t+1 (Ft+2 (x)) .
Therefore, the variance of Fb(x) satisfies TX −1 2 21 1 t+1 t T T V Fb(x) 2 ≤E Fb1 (Ft+1 (x)) − Fb1 (Ft+2 (x)) ≤ =
t=0 T −1 X
E
t=0 T −1 X t=0
2 21 t+1 T t T b b F1 (Ft+1 (x)) − F1 (Ft+2 (x))
1 2 2 t t T T b t ft+1 ξt+1 , F (x) , E Fb1 Et ft+1 ξt+1 , Ft+2 (x) − Fb1 E t+2
where the second inequality follows from the Minkowski inequality. From now on, we use the notational T (x)) for all t ∈ [T ]. As F bt is L[t] -Lipschitz continuous, we thus find shorthand At := ft (ξt , Ft+1 1 h −1 i2 12 1 TX 2 b b L[t] E E Et [At+1 ] − Et [At+1 ] V Fb(x) ≤ ≤ ≤
t=0 T −1 X t=0 T −1 X t=0
2
1 2 2 b L[t] E Et [At+1 ] − Et [At+1 ] 2
1 TX −1 2 2 L[t] σt+1 L[t] ≤ E A − E [A ] . √ √ t+1 t t+1 nt+1 nt+1 2 t=0
b Here, the second inequality exploits Jensen’s inequality and the observation that E[E[·]] = E[·]. The third inequality follows from the formula for the (conditional) variance of a sample average, and the last inequality follows from Assumption 3.5. Squaring both sides of the above inequality then yields the basic bound TX −1 L[t] σt+1 2 b . (29) V F (x) ≤ √ nt+1 t=0
P This basic bound can be strengthened by recalling that Fb(x) = n11 ni11=1 Fbi1 (x), where Fbi1 (x), i1 ∈ [n1 ], constitute independent and identically distributed random variables that depend on the samples of the i1 -th scenario tree used in the construction of the SAA estimator; see Definition 3.4. Thus, we have T −1 2 n1 1 X L[t] σt+1 1 X b b V Fi1 (x) ≤ + σ1 V F (x) = 2 √ n1 nt+1 n1 t=1
i1 =1
where the inequality uses (29) with n1 = 1 to bound the variance of Fbi1 (x). Hence, the claim follows. 42
The next lemma shows that sub-Gaussianity is preserved under averaging and concatenation as well as under transformations with Lipschitz continuous functions and mixtures. Among other things, this lemma enables us to identify easily verifiable conditions under which Assumption 3.11 holds. Lemma C.5 (Operations Preserving Sub-Gaussianity). The following hold. (i) If zi , i ∈ P [n], are i.i.d. sub-Gaussian random vectors with common variance proxy ζ 2 , then their 1 average n ni=1 zi constitutes another sub-Gaussian random vector with variance proxy ζ 2 /n. (ii) If z1 and z2 are (possibly dependent) sub-Gaussian random vectors valued in Rm1 and Rm2 with respective variance proxies ζ12 and ζ22 , then (z1 , z2 ) is sub-Gaussian with variance proxy 2 max{ζ12 , ζ22 }. (iii) If z is an m-dimensional sub-Gaussian random vector with variance proxy ζ 2 , and if f : Rm → Rn is an L-Lipschitz continuous function with ∥f (z) − f (z ′ )∥2 ≤ L∥z − z ′ ∥2 for all z, z ′ ∈ Rm , then f (z) is a sub-Gaussian random vector with variance proxy that grows monotonically with L2 m2 ζ 2 . (iv) If zi , i ∈ [n] are (possibly dependent) sub-Gaussian random vectors valued in Rm with respective variance proxies ζi2 , i ∈ [n], and if α is a random index independent of all other random objects with P(α = i) = wi for all i ∈ [n], then the mixture zα is a sub-Gaussian random vector with variance proxy maxi∈[n] {ζi2 } + maxi∈[n] {∥E[zi ] − E[zα ]∥2 }.
Proof. Assertion (i) follows from an immediate generalization of [64, Exercise 2.13]. As for assertion (ii), observe that ⊤ E exp (λ1 , λ2 )⊤ ((z1 , z2 ) − E[(z1 , z2 )]) = E exp λ⊤ 1 (z1 − E[z1 ]) exp λ2 (z2 − E[z2 ]) 1/2 1/2 ≤ E exp 2λ⊤ E exp 2λ⊤ 1 (z1 − E[z1 ]) 2 (z2 − E[z2 ])
≤ exp(∥λ1 ∥22 ζ12 ) exp(∥λ2 ∥22 ζ22 ) ≤ exp(∥(λ1 , λ2 )∥22 max{ζ12 , ζ22 })
for all λ1 ∈ Rm1 and λ2 ∈ Rm2 . Here, the three inequalities follow from Hölder’s inequality, the subGaussianity of z1 and z2 and the monotonicity of the exponential function, respectively. This shows that the combined random vector (z1 , z2 ) is indeed sub-Gaussian with variance proxy 2 max{ζ12 , ζ22 }. As for assertion (iii), let z ′ be an independent copy of z. Thus, −z ′ has the same variance proxy as z, which implies via [64, Exercise 2.13] that z − z ′ is sub-Gaussian with variance proxy 2ζ 2 . Hence, we find E exp λ⊤ (f (z) − E[f (z)]) = E exp λ⊤ (E[f (z) − f (z ′ )|z]) ≤ E exp λ⊤ f (z) − f (z ′ ) ≤ E exp ∥λ∥2 ∥f (z) − f (z ′ )∥2 ≤ E exp L ∥λ∥2 ∥z − z ′ ∥2 , where we have used Jensen’s inequality, the Cauchy-Schwarz inequality and the L-Lipschitz continuity of f . Since ∥ · ∥2 ≤ ∥ · ∥1 and as the exponential function is monotonically increasing, we also find exp(∥z − z ′ ∥2 ) ≤ exp(∥z − z ′ ∥1 ) =
m Y i=1
exp(|zi − zi′ |).
Combining the above estimates and applying the multi-variable Hölder inequality then yields Y m ⊤ ′ E exp λ (f (z) − E[f (z)]) ≤ E exp L ∥λ∥2 |zi − zi | i=1
m Y 1/m ≤ E exp Lm ∥λ∥2 |zi − zi′ | . i=1
43
(30)
As z 1/m increases with z ≥ 0, an upper bound on (30) follows by bounding the underlying expectations as E exp Lm ∥λ∥2 |zi − zi′ | ≤ E exp Lm ∥λ∥2 (zi − zi′ ) + E exp Lm ∥λ∥2 (zi′ − zi ) ≤ 2 exp ∥λ∥22 L2 m2 ζ 2 . Here, the two inequalities hold because exp(|a|) ≤ exp(a) + exp(−a) for any a ∈ R and because zi − zi′ is sub-Gaussian with variance proxy 2ζ 2 . Substituting this estimate into (30) shows that E exp λ⊤ (f (z) − E[f (z)]) ≤ 2 exp ∥λ∥22 L2 m2 ζ 2 . (31) It remains to be shown that the factor 2 can be absorbed in the variance proxy. For ease of notation, we henceforth use ψ(λ) as shorthand for the left hand side of (31). That is, ψ(λ) denotes the moment-generating function of f (z) − E[f (z)]. We must now prove that there exists ς > 0 such that ψ(λ) ≤ exp(∥λ∥22 ς 2 /2) for all λ ∈ Rn . To this end, note that (31) ensures that ψ(λ) is analytic on a neighborhood of 0 [64, Section 2]. Thus, there exists ς12 > 0 with ∇2 ψ(0) ⪯ ς12 Im . In addition, we have ψ(0) = 1 and ∇ψ(0) = 0. Next, define ϕ1 (λ) := exp(∥λ∥2 ς12 /2), and note that ϕ1 (0) = 1, ∇ϕ1 (0) = 0 and ∇2 ϕ1 (0) = ς12 Im . Thus, Taylor’s theorem guarantees that there exists ϵ > 0 such that ψ(λ) ≤ ϕ1 (λ) whenever ∥λ∥2 ≤ ϵ. Define now ς22 := 2(log(2)/ϵ2 + L2 m2 ζ 2 ), and set ϕ2 (λ) := exp(∥λ∥2 ς22 /2). By (31), we then have ψ(λ) ≤ exp log(2) + ∥λ∥22 L2 m2 ζ 2 = exp ∥λ∥22 log(2)/∥λ∥22 + L2 m2 ζ 2 < ϕ2 (λ) whenever ∥λ∥2 > ϵ. We may thus conclude that ψ(λ) ≤ max{ϕ1 (λ), ϕ2 (λ)} = exp ∥λ∥2 max{ς12 , ς22 }/2
∀λ ∈ Rm .
Consequently, f (z) is sub-Gaussian with variance proxy max{ς12 , ς22 }. As for assertion (iv), observe that n X ⊤ E exp(λ (zα − E[zα ])) = wi E exp(λ⊤ (zi − E[zα ])) i=1
=
n X i=1
wi exp(λ⊤ (E[zi ] − E[zα ])) E exp(λ⊤ (zi − E[zi ])) ,
where the first equality holds because α is independent of zi , i ∈ [n], and the second equality centers each zi around its own mean. The sub-Gaussianity of each zi implies that E[exp(λ⊤ (zi −E[zi ]))] ≤ exp(∥λ∥22 ζi2 /2). 2 Factoring out a uniform bound involving ζmax := maxi∈[n] {ζi2 }, we obtain n X ⊤ 2 2 E exp(λ (zα − E[zα ])) ≤ exp(∥λ∥2 ζmax /2) wi exp(λ⊤ (E[zi ] − E[zα ])). i=1
Here, the sum can be interpreted as the moment generating function of aP discrete random variable Y taking values yi = λ⊤ (E[zi ] − E[zα ]) with probabilities wi . Note that E[Y ] = wi λ⊤ (E[zi ] − E[zα ]) = 0. Note also that |yi | ≤ ∥λ∥2 R for R := maxi∈[n] {∥E[zi ]−E[zα ]∥2 }, that is, the variable Y is bounded. Hoeffding’s Lemma (see, e.g., [64, Example 2.4]) thus implies that E[exp(Y )] ≤ exp(∥λ∥22 R2 /2). In summary, we find 2 E exp(λ⊤ (zα − E[zα ])) ≤ exp ∥λ∥22 (ζmax + R2 )/2 , 2 which establishes that zα is sub-Gaussian with variance proxy ζmax + R2 .
44
Proof of Lemma 3.13. Assumptions 3.1 and 3.2 ensure that X is compact and that both F and Fb are L[T ] Lipschitz continuous; see Lemma C.3 as well as the discussion after Definition C.4. Next, set v = ϵ/(4L[T ] ). By Lemma A.1, there exists a v-net Xv ⊆ X with |Xv | ≤ ⌈2DX /v + 1⌉d such that X ϵ . (32) P sup |Fb(x) − F (x)| > ϵ ≤ P |Fb(x′ ) − F (x′ )| > 2 x∈X ′ x ∈Xv
P Recall that the SAA estimator is given by Fb (x) = n11 ni11=1 Fbi1 (x), where the Fbi1 (x), i1 ∈ [n1 ], can be interpreted as independent and identically distributed estimators for F (x), each of which is constructed from a single scenario tree; see Definition 3.4. In the following we introduce the auxiliary bias function µ(x) := E Fbi1 (x) − F (x) = E Fb(x) − F (x) . Given any ϵ > 0, we now seek sufficiently large sample sizes to ensure that |µ(x′ )| ≤ 4ϵ for every x′ ∈ Xv . If the integrands are Lipschitz continuous but not smooth (i.e., only Assumption 3.2 holds), then we have 4L[t] σt+1 (T − 1) 2 L[t] σt+1 ϵ nt+1 = ∀t ∈ [T − 1] =⇒ ∀t ∈ [T − 1], ≤ √ ϵ nt+1 4(T − 1)
which ensures via Lemma 3.6 that
µ(x′ ) ≤
T −1 X t=1
L[t] σt+1 ϵ ≤ √ nt+1 4
∀x′ ∈ Xv .
Similarly, if the integrands are Lipschitz continuous and smooth (i.e., both Assumptions 3.2 and 3.3 hold), then we have 2 (T − 1) 2 2L[t−1] St σt+1 L[t−1] St σt+1 ϵ nt+1 = ∀t ∈ [T − 1] =⇒ ≤ ∀t ∈ [T − 1], ϵ 2nt+1 4(T − 1) which ensures via Lemma 3.6 that
µ(x′ ) ≤
T −1 X t=1
2 L[t−1] St σt+1 ϵ ≤ 2nt+1 4
∀x′ ∈ Xv .
In both cases, we may conclude that ϵ ϵ P Fb(x′ ) − F (x′ ) > ≤ P Fb(x′ ) − F (x′ ) > + µ(x′ ) 2 4 ! n 1 n1 ϵ2 1 X b ϵ ′ ′ ′ =P Fi1 (x ) − F (x ) − µ(x ) > ≤ exp − n1 4 32ζ 2 i1 =1
for all x′ ∈ Xv . The second inequality follows from a standard Hoeffding bound [64, Proposition 2.5], which applies because the i.i.d. centered random variables Fbi1 (x) − F (x) − µ(x), i ∈ [n1 ], are sub-Gaussian with variance proxy ζ 2 thanks to Assumption 3.11. Following a similar reasoning, one can also show that ϵ n1 ϵ2 ′ ′ b P F (x ) − F (x ) > ≤ exp − 2 32ζ 2 for all x′ ∈ Xv . Combining the above probability bounds with (32) finally yields d n1 ϵ 2 n1 ϵ2 b P sup |F (x) − F (x)| > ϵ ≤ 2|Xv | exp − ≤ 2 8L[T ] DX /ϵ + 1 exp − , 32ζ 2 32ζ 2 x∈X and thus the claim follows. 45
Proof of Corollary 3.15. Problems (MCCO) and (1) are solvable by virtue of Weierstrass’ extreme value theorem, which applies thanks to Assumptions 3.1 and 3.2 and Lemma C.3. Hence, x∗ and x b∗ exist. Elementary manipulations then show that P (F (b x∗ ) − F (x∗ ) > ϵ) = P F (b x∗ ) − Fb(b x∗ ) + Fb(b x∗ ) − Fb(x∗ ) + Fb(x∗ ) − F (x∗ ) > ϵ ≤ P F (b x∗ ) − Fb(b x∗ ) > ϵ/2 or Fb(x∗ ) − F (x∗ ) > ϵ/2 ≤ P F (b x∗ ) − Fb(b x∗ ) > ϵ/2 + P Fb(x∗ ) − F (x∗ ) > ϵ/2 d n1 ϵ2 ≤ 4 8L[T ] DX /ϵ + 1 exp − , 128ζ 2 where the first inequality holds because Fb(b x∗ ) − Fb(x∗ ) ≤ 0. The second and the third inequalities follow from the union bound and from Lemma 3.13, respectively. Proof of Theorem 3.16. The condition on n1 ensures via Corollary 3.15 that P (F (b x∗ ) − F (x∗ ) > ϵ) ≤ ∗ ∗ β for any minimizers x and x b of problems (MCCO) and (1), respectively. Hence, any minimizer of problem (1) is an ϵ-optimal solution to (MCCO) with probability at least 1 − β. The scenario complexity of the SAA estimator is thus obtained by multiplying the given expressions for the sample sizes nt , t ∈ [T ]. Next, we provide the proofs of all formal results in Section 3.2. Proof of Lemma 3.19. Since x is fixed, we simplify notation by suppressing the dependence on x for all functions of x throughout the proof. Suppose first that only Assumptions 3.1, 3.2 and 3.18 hold. We prove the upper bound on µ2t by backward induction on t. As for the base case corresponding to t = T , note that b T ∥22 ] = ET −1 [∥fT (ξT , x)∥22 ] ≤ µ2T , µ2T = ET −1 [∥H where the inequality follows from the definition of µ2T . As for the induction step corresponding to any t b t (x) in (18) to obtain with t + 1 ∈ [T ], we use the law of total expectation and the definition of H b t ∥2 µ2t = Et−1 ∥H 2
=
=
Mt X ℓ=0 Mt X ℓ=0
qt (ℓ) Et−1
1 1 b λt ,e 1 b λt ,o 2 λt b h − ht − ht λt = ℓ 2 qt (λt )2 t 2 2
(33)
h 1 1 ℓ,e 1 b ℓ,o 2 i Et−1 b h − ht 2 . hℓt − b qt (ℓ) 2 t 2
In the second line we have eliminated the condition on the event λt = ℓ, which is allowed because λt is 1 b ℓ,o independent of all other random objects. Next, define ∆ℓ := b hℓt − 12 b hℓ,e t − 2 ht to avoid clutter, and note that Lt b ℓ,e b b ℓ [H b t+1 ] + Et [Ht+1 ] − E t 2 2 Lt b ℓ,e b b ℓ,o [H b t+1 ] = Et [Ht+1 ] − E t 2 2 Lt b ℓ,e b b t+1 ] + ≤ Et [Ht+1 ] − Et [H 2 2
∥∆ℓ ∥2 ≤
b ℓ,o [H b ℓ [H b t+1 ] − E b t+1 ] E t t 2
b ℓ,o [H b t+1 ] − Et [H b t+1 ] , E t 2
where the first inequality holds because of Assumption 3.2, which ensures that ft (ξt , xt ) is Lt -Lipschitz continuous in xt , and because of the triangle inequality. The equality follows from the antithetic sample
46
average identity (17), and the second inequality follows again from the triangle inequality. As (a + b)2 ≤ 2a2 + 2b2 for all a, b ∈ R, the conditional second moment of ∆ℓ admits the following upper bound. h i h i L2 b ℓ,e [H b ℓ,o [H b t+1 ] − Et [H b t+1 ] 2 + Et−1 E b t+1 ] − Et [H b t+1 ] 2 Et−1 ∥∆ℓ ∥22 ≤ t Et−1 E t t 2 2 2 i B L2 h 2 b ℓ,e [H b t+1 ] − Et [H b t+1 ] 2 ≤ 2 t Et−1 ∥H b t+1 ∥22 = B2 Lt Et−1 µ2t+1 . = L2t Et−1 E t 2 ℓ−1 ℓ−1 2 2 b t have the same distribution Here, the first equality holds because the even and the odd sample averages of H conditional on Ft , and the second inequality follows from Lemma A.4. The second equality follows from the definition of µ2t+1 and the law of total expectation. Substituting this bound back into (33) yields µ2t ≤ 2B2 L2t Et−1 µ2t+1
≤ 2B2 L2t µ2T
TY −1 s=t+1
Mt X ℓ=0
1 qt (ℓ)2ℓ
2B2 L2s
Ms X ℓ=0
1 s 2 qs (ℓ)
X Mt ℓ=0
Ms TY −1 X 1 1 2 2 = µT 2B2 Ls , 2s qs (ℓ) qt (ℓ)2ℓ s=t ℓ=0
where the second inequality follows from the induction hypothesis. This establishes the first claim. t Suppose now that Assumption 3.3 holds, too. In this case, we prove the upper bound on µ2t by backward induction on t. The base case corresponding to t = T is trivially true because T b T ∥2TT = ET −1 ∥fT (ξT , x)∥2TT ≤ µ2T , µ2T = ET −1 ∥H T 2 2 T
where the inequality uses the definition of µ2T . As for the induction step corresponding to any t with t + 1 ∈ [T ], we can use the law of total expectation and define ∆ℓ as in the first part of the proof to demonstrate that t b t ∥2tt ] = µ2t = Et−1 [∥H 2
t Mt X Et−1 ∥∆ℓ ∥2t 2
ℓ=0
qt (ℓ)2t −1
.
(34)
Next, we can use the antithetic sample average identity (17) to reformulate ∆ℓ as b ℓ [H b t+1 ] − Et [H b t+1 ] b t+1 ]) − ∇xt ft (ξt , Et [H b t+1 ])⊤ E ∆ℓ = b hℓt − ft (ξt , Et [H t i 1 h ℓ,e b ℓ,e [H b t+1 ]) − ∇xt ft (ξt , Et [H b t+1 ])⊤ E b t+1 ] − Et [H b t+1 ] − b ht − ft (ξt , Et [H t 2 h i 1 ℓ,o b ℓ,o [H b t+1 ])⊤ E b t+1 ] − Et [H b t+1 ] . b t+1 ]) − ∇xt ft (ξt , Et [H − b ht − ft (ξt , Et [H t 2 b ℓ,o Taking norms on both sides and replacing b hℓt , b hℓ,e t and ht by their definitions then yields ∆ℓ 2t ≤ ∆ℓ 2 ≤
2 St b ℓ b b t+1 ] Et [Ht+1 ] − Et [H 2 2 2 2 St b ℓ,e b b ℓ,o [H b t+1 ] + St E b t+1 ] − Et [H b t+1 ] . + Et [Ht+1 ] − Et [H t 4 4 2 2
Here, the two inequalities follow from Lemma A.2 and from a combination of the triangle inequality and
47
the inequality (13) for smooth functions, respectively. By the Minkowski inequality, we then obtain t 1 Et−1 ∥∆ℓ ∥22t 2t h i 1t h i 1t St St 2t+1 2t+1 2 2 ℓ,e b ℓ b b b b b ≤ Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 + Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 2 4 i 1t h St 2t+1 2 ℓ,o b b b + Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 4 h h i 1t i 1t St St 2t+1 2t+1 2 2 ℓ,e b ℓ b b b b b + , = Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 2 2 b t have the same distribution where the equality holds because the even and the odd sample averages of H t+1 conditional on Ft . Applying Lemma A.2 with p = 2, q = 2 , m = 1 and n = dt to the two differences of conditional expectations in the above expression then yields t 1 Et−1 ∥∆ℓ ∥22t 2t i 1t h St 2t+1 2t −1 ℓ b b b (dt ) Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2t+1 2 ≤ 2 h i 1t St 2t+1 ℓ,e b 2t −1 b b (dt ) Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2t+1 2 + 2 h h i 1t i 1t St St 2t+1 2t+1 2 2t −1 B2t+1 2t −1 B2t+1 b b (dt ) + (dt ) ≤ Ht+1 2t+1 2 Ht+1 2t+1 t ℓ Et−1 t (ℓ−1) Et−1 2 2 2 2 2 2 h i 1t h t+1 i 1t t+1 3St B 3S t t t+1 2 t 2 2 2 −1 B2t+1 b t+1 t+1 = (dt )2 −1 2t ℓ Et−1 H = (d ) E µ2t+1 2 . t t−1 t 2 2 ℓ 2 2 2 2 Here, the second inequality follows from Lemma A.4 with p = 2t+1 and n = 2ℓ (for sample averages involving even as well as odd samples) or n = 2ℓ−1 (for sample averages involving only even or odd samples). Next, by raising both sides of the resulting inequality to the 2t -th power, we find t h t+1 i ℓ 2t 3St 2 1 t Et−1 ∥∆ ∥2t ≤ (dt )2 −1 B2t+1 2t ℓ Et−1 µ2t+1 2 2 2t s Ms TY −1 X 1 1 3St 3Ss 2 2s −1 t T (dt )2 −1 B2t+1 2t ℓ µ2T ds B2s+1 ≤ , 2s ℓ′ q (ℓ′ )2s −1 2 2 2 2 s ′ s=t+1 ℓ =0
where the second inequality uses the induction hypothesis. Substituting this back into (34) finally yields t
T
µ2t ≤ µ2T
TY −1 s=t
3Ss 2
2s
s
d2s −1 B2s+1
Ms X ℓ=0
1 . 22s ℓ qs (ℓ)2s −1
This observation completes the proof. b t (x) represents an estimator for the conditional risk mapping F T (x) introduced The random variable H t b t (x) for any t ∈ [T ]. in Definition C.1. Lemma C.6 below bounds the conditional bias of the estimator H b 1 (x)] for any x ∈ X , Lemma 3.20 follows immediately from Lemma C.6. Its proof is As E[Fb(x)] = E[H thus omitted for brevity. We point out that the stronger result of Lemma C.6 will also be used in Section 4. Lemma C.6. If Assumptions 3.1, 3.2 and 3.18 hold and if t ∈ [T ], then we have 1 PT −1 L[t:s] Et−1 [µ2s+1 (x)] 2 if Assumption 3.3 does not hold, T s=t b t (x)] − Ft (x) ≤ 2Ms /2 Et−1 [H 2 PT −1 L[t:s−1] Ss Et−1 [µ2s+1 (x)] if Assumption 3.3 holds. s=t 2Ms +1 48
Proof. We simplify notation by suppressing the dependence on x for all functions of x throughout the proof. Suppose first that only Assumptions 3.1 and 3.2 hold. By the construction of the conditional risk mappings b s in Definition 3.17, we have Fst in Definition C.1 and the estimators H b T ] and Et−1 [H b t ] = F t−1 Et−1 [H bt] . FtT = FtT −1 ET −1 [H t b t admits the estimate The triangle inequality thus implies that the squared conditional bias of H b t ] − FtT Et−1 [H
≤ 2
T −1 X s=t
b s ] − Fts Es [H b s+1 ] . Fts−1 Es−1 [H 2
To reformulate the resulting upper bound, observe that for any s ∈ [T − 1], we have 1 1 b λs ,e 1 b λs ,o λs b b Es−1 [Hs ] = Es−1 hs − hs − hs qs (λs ) 2 2 M s X 1 ℓ,e 1 b ℓ,o s b Ms b hs − hs = Es−1 [b hM = Es−1 b hℓs − b s ] = Es−1 fs (ξs , Es [Hs+1 ]) , 2 2
(35)
(36)
ℓ=0
b s in (18), and the second uses the law of total where the first equality follows from the definition of H expectation along with the independence of λs from all other random objects. The third and fourth equalb ℓ,o b ℓ−1 conditional on Fs for all ℓ ∈ [Ms ] and ities hold because b hℓ,e s and hs share the distribution of hs 0,e 0,o b b because hs = hs = 0, respectively. Hence, the s-th term in the sum on the right hand side of (35) satisfies b s ] − F s Es [H b s+1 ] Fts−1 Es−1 [H t 2 b Ms [H b b s+1 ])] = Fts−1 Es−1 [fs (ξs , E ])] − Fts−1 Es−1 [fs (ξs , Es [H s+1 s 2 Ms b b b ≤ L[t:s] Et−1 Es [Hs+1 ] − Es [Hs+1 ] 2 Ms 1 b [H b s+1 ] − Es [H b s+1 ] 2 2 ≤ L[t:s] Et−1 E s 2 1 1 L[t:s] L 2 b s+1 − Es [H b s+1 ] 2 2 ≤ [t:s] Et−1 µ2 ≤ M /2 Et−1 H , s+1 2 2 s 2Ms /2 where the first inequality follows from Lemma C.3, which guarantees that Fts−1 (zs−1 ) is L[t:s−1] -Lipschitz continuous in zs−1 , and from Assumption 3.2, which ensures that fs (ξs , xs ) is Ls -Lipschitz continuous in xs . The second and third inequalities follow from Jensen’s inequality and the fact that the variance of the average of i.i.d. random variables equals the variance of one of these random variables divided by the sample size. The last inequality holds because the variance of any random variable is bounded by its second moment. The claim then follows by substituting the emerging bound into (35). If Assumption 3.3 holds, too, then the s-th term on the right hand side of (35) satisfies b s ] − Fts Es [H b s+1 ] Fts−1 Es−1 [H 2 s−1 Ms b b b s+1 ])] (37) = Ft Es−1 [fs (ξs , Es [Hs+1 ])] − Fts−1 Es−1 [fs (ξs , Es [H 2 b Ms [H b s+1 ]) − fs (ξs , Es [H b s+1 ]) . ≤ L[t:s−1] Et−1 Es−1 fs (ξs , E s 2 b Ms [·]] = Es [·], we may conclude that b s+1 ]) is Fs -measurable and as Es [E As ∇xs fs (ξs , Es [H s b Ms [H b s+1 ])⊤ (E b s+1 ] − Es [H b s+1 ]) = 0. Es−1 ∇xs fs (ξs , Es [H s 49
Hence, the last expression in (37) equals b Ms [H b s+1 ]) − fs (ξs , Es [H b s+1 ]) L[t:s−1] Et−1 Es−1 fs (ξs , E s b Ms [H b s+1 ])⊤ (E b s+1 ] − Es [H b s+1 ]) − ∇xs fs (ξs , Es [H s
2
Ms L[t:s−1] Ss b [H b s+1 ] − Es [H b s+1 ] 2 ≤ Et−1 E s 2 2 L S L[t:s−1] Ss b s+1 − Es [H b s+1 ] 2 ≤ [t:s−1] s Et−1 µ2s+1 . ≤ Et−1 H 2 M +1 M +1 s s 2 2
Here, the first inequality follows by applying Jensen’s inequality to move the conditional expectation Es−1 [·] outside the norm, and then using the smoothness inequality (13), which holds by Assumption 3.3. The second inequality holds because the variance of an average of i.i.d. samples equals the variance of a single sample divided by the sample size. The third inequality follows from the fact that the variance of any random variable is bounded by its second moment, together with the definition of µ2s+1 in (19). Substituting the emerging bound into (35) yields the desired result. Proof of Lemma 3.22. The MLMC estimator of Definition 3.17 is built on n1 independent scenario trees, b 1 (x))]. Each tree accommodates T stages, 2λt branches per and thus, we have E[C(Fb(x))] = n1 E[C(H stage-t node for every t ∈ [T − 1] and one sample per node. As each log-branching factor λt , t ∈ [T − 1], is independent of all other random objects, the expected number of scenarios per tree is given by Mt −1 X TY b 1 (x)) = E C(H qt (ℓ)2ℓ . t=1 ℓ=0
Hence, the claim follows. The following lemma shows that if the integrands fail to be smooth, then the rate parameters in Assumption 3.23 are optimal in a precise sense, as detailed in Remark 3.24. Lemma C.7. If Assumptions 3.1, 3.2 and 3.18 hold, then the choice r = (r1 , . . . , rT −1 ) with rt = 1/2 for all t ∈ [T − 1] minimizes the approximate sampling cost cns (r) over all r ∈ (0, 1)T −1 . Proof of Lemma C.7. We adopt a similar proof strategy as in [11, § 3], which proposes an untruncated MLMC estimator for a subclass of (MCCE) problems with T = 2. Using the definition of cns (r) and the uniform upper bound on µ2t (x) for nonsmooth integrands established in Lemma 3.19, we obtain ! T −1 M ! ! M ! Mt T −1 Mt TY −1 X t t YX X 2 YX 1 2 1 cns (r) = 2 qt (ℓ)2ℓ = 2 qt (ℓ)2ℓ . ϵ ϵ qt (ℓ)2ℓ qt (ℓ)2ℓ t=1 ℓ=0
t=1 ℓ=0
t=1
ℓ=0
ℓ=0
As qt (ℓ) = rt (1 − rt )ℓ /(1 − (1 − rt )Mt +1 ) is independent of rt′ for all t′ ̸= t, we can find the optimal rt by minimizing the logarithm of the t-th term in the resulting product over rt ∈ (0, 1). Thus, we must solve min log
rt ∈(0,1)
X Mt ℓ=0
1 ℓ 2 qt (ℓ)
+ log
X Mt
2 qt (ℓ) , ℓ
ℓ=0
where qt (ℓ) is interpreted as a function of rt . This minimization problem can be simplified by representing the rate as rt = 1 − 2−bt using an auxiliary decision variable bt ∈ (0, ∞). Under this re-parametrization
50
we have qt (ℓ) = 2−bt ℓ /Zt , and the normalization constant amounts to Zt = (1 − 2−bt (Mt +1) )/(1 − 2−bt ). Replacing qt (ℓ) with this formula allows us to reformulate the above minimization problem equivalently as min
bt ∈(0,∞)
log
X Mt
(bt −1)ℓ
2
+ log
ℓ=0
X Mt
2
(1−bt )ℓ
.
ℓ=0
The objective function of the re-parametrized optimization problem is convex in bt because logarithms of sums of exponential functions are known to be convex. An optimal bt can thus be found from the problem’s first-order optimality condition. Indeed, the gradient of the objective function vanishes if and only if PMt PMt (bt −1)ℓ ℓ 2(1−bt )ℓ ℓ=0 ℓ 2 . PMt (b −1)ℓ′ = Pℓ=0 Mt (1−bt )ℓ′ t ℓ′ =0 2 ℓ′ =0 2 Clearly, the above equation is solved by bt = 1 for every t ∈ [T − 1], which implies that cns (r) is minimized by setting rt = 1 − 2−bt = 1/2 for every t ∈ [T − 1]. Hence, the claim follows. Proof of Theorem 3.25. To guarantee a root mean squared error of at√most ϵ, we require both the bias and the standard deviation of the MLMC estimator to be bounded by ϵ/ 2 for any x ∈ X . Suppose first that Assumption 3.3 does not hold (i.e., the integrands are only Lipschitz continuous). Hence, the bias satisfies E[Fb(x)] − F (x) ≤
1
T −1 X
L[t] E[µ2t+1 (x)] 2
t=1
2Mt /2
thanks to Lemma 3.20. Lemma 3.19 further implies that 1 E[µ2t+1 (x)] 2 ≤
p
Ct+1
Ms TY −1 X s=t+1
ℓ=0
1 ℓ 2 qs (ℓ)
1 2
=
TY −1 p p p Ct+1 Z s Ms + 1 s=t+1
∀t ∈ [T − 1],
where the equality holds thanks to Assumption 3.23, which implies that rt = 1/2. Indeed, for this choice of the rate parameter rt we have qt (ℓ) = 2−ℓ /Zt , and the corresponding normalization constant simplifies to Zt = 2(1 − 2−(Mt +1 ). Combining the above estimates then yields E[Fb(x)] − F (x) ≤
√ −1 p p L[t] Ct+1 TY Zs Ms + 1. 2Mt /2 s=t+1 t=1
T −1 X
Recalling from the theorem statement that the t-th truncation point is given by √ Q √ T −1 √ √ Mt = 2 log2 2L[t] Ct+1 s=t+1 Zs Ms + 1 (T − 1)/ϵ √ for every t ∈ [T − 1], we finally obtain the desired estimate |E[Fb(x)] − F (x)| ≤ ϵ/ 2 for every x ∈ X . Similarly, the variance of the MLMC estimator satisfies 2 2 bH b 1 (x)] = 1 V H b 1 (x) ≤ µ1 (x) ≤ ϵ V Fb(x) = V E[ n1 n1 2
b respectively. The first for every x ∈ X , where the two equalities exploit the definitions of Fb(x) and E[·], inequality holds because the variance of any random variable is bounded by its second moment, and the second inequality follows from our choice of n1 = ⌈supx∈X 2µ21 /ϵ2 ⌉. In summary, the above reasoning implies that the root mean squared error of Fb(x) is bounded by ϵ uniformly across all x ∈ X . 51
Next, Lemma 3.22 implies that the expected sampling cost of the MLMC estimator amounts to E[C(Fb(x))] = n1
Mt TY −1 X
qt (ℓ)2ℓ = n1
TY −1
t=1 ℓ=0
Zt−1 (Mt + 1),
t=1
where the second equality follows from Assumption 3.23. It remains to analyze how the resulting expression scales with respect to 1/ϵ. The recursive definition of the truncation points readily implies that Mt = O(log(ϵ−1 )) for all t ∈ [T − 1]. By Lemma 3.19 and the definition of n1 , we further have Ms T −1 X
2C1 Y n1 ≤ 2 ϵ
s=1
ℓ=0
1 ℓ 2 qs (ℓ)
T −1
2C1 Y +1= 2 Zt (Mt + 1) + 1 = O(log(ϵ−1 )T −1 ϵ−2 ). ϵ t=1
Hence, the expected sampling cost of the MLMC estimator satisfies E[C(Fb(x))] ≤ O log(ϵ−1 )2(T −1) ϵ−2 . Suppose now that Assumption 3.3 holds, too. In this case, Lemma 3.20 and Jensen’s inequality imply that the bias of Fb(x) satisfies E[Fb(x)] − F (x) ≤
T −1 X t=1
1
−1 L[t−1] St E[µ2t+1 (x)] 2t L[t−1] St E[µ2t+1 (x)] TX ≤ . 2Mt +1 2Mt +1 t+1
t=1
Lemma 3.19 and Assumption 3.23 further imply that 1 t+1 E[µ2t+1 (x)] 2t ≤
Dt+1
Ms TY −1 X
ℓ s Zs2 −1 2− 2s
1t 2
TY −1
2−t = Dt+1
2s −1 2t
Zs
1−2
s=t+1
s=t+1 ℓ=0
Ms +1
1 − 2 − 2s
1t 2
− 21s
−t
∀t ∈ [T − 1].
−t
Indeed, setting the t-th rate parameter to rt = 1 − 2−1−2 implies that qt (ℓ) = 2−(1+2 )ℓ /Zt with normal−t −t ization constant Zt = (1 − 2−(1+2 )(Mt +1) )/(1 − 2−1−2 ). Combining these estimates then yields E[Fb(x)] − F (x) ≤
T −1 X t=1
−t
TY −1 2 2s −1 L[t−1] St Dt+1 2t Z s 2Mt +1
s=t+1
Ms +1
1 − 2 − 2s 1−2
1t 2
− 21s
.
Recall now from the theorem statement that the t-th truncation point is given by 1t 2s −1 √ − Mss+1 QT −1 −t 2 2 1−2 2 2t Mt = log2 2L[t−1] St Dt+1 (T − 1)/(2ϵ) s=t+1 Zs − 1 1−2
2s
√ for every t ∈ [T − 1]. This choice then yields the desired bias estimate |E[Fb(x)] − F (x)| ≤ ϵ/ 2 for every x ∈ X . To bound the variance of the MLMC estimator, we may use the exact same reasoning as in the first part of the proof. Thus, we may conclude that V(Fb(x)) ≤ ϵ2 /2 uniformly across all x ∈ X . In addition, we may use a similar reasoning as in the first part of the proof to conclude that E[C(Fb(x))] = n1
Mt TY −1 X
−t Zt−1 2−2 ℓ = n1
t=1 ℓ=0
TY −1 t=1
M +1
− 2tt −1 1 − 2 Zt − 21t
1−2
.
The recursive definition of the truncation points readily implies that Mt = O(log(ϵ−1 )) for all t ∈ [T − 1]. Furthermore, the definition of n1 implies via Lemma 3.19 that T −1 M
T −1
Mt +1
t t 2D1 Y X 2D1 Y 2t −1 1 − 2− 2t −t n1 ≤ 2 Zt2 −1 2−2 ℓ + 1 = 2 Zt 1 ϵ ϵ 1 − 2 − 2t
t=1 ℓ=0
t=1
+ 1 = O(ϵ−2 ).
Hence, the expected sampling cost of the MLMC estimator satisfies E[C(Fb(x))] ≤ O(ϵ−2 ). 52
Proof of Theorem 3.31. The condition on n1 ensures via Corollary 3.30 that P (F (b x∗ ) − F (x∗ ) > ϵ) ≤ β for any minimizers x∗ and x b∗ of problems (MCCO) and (1), respectively. Hence, any minimizer of problem (1) constitutes an ϵ-optimal solution to (MCCO) with probability at least 1 − β. By Lemma 3.22, the expected sampling cost of the MLMC estimator amounts to E[C(Fb(x))] = n1
Mt TY −1 X
qt (ℓ)2ℓ .
t=1 ℓ=0
By using a similar reasoning as in the proof of Theorem 3.25, one can show that the product term is of order O(log(ϵ−1 )T −1 ) for Lipschitz continuous integrands and of order O(1) for smooth integrands. Moreover, the number n1 of scenario trees appearing in the statement of the corollary is of order O(log(ϵ−1 )ϵ−2 ). The desired scenario complexity bounds then follow by multiplying these estimates.
D
Proofs of Section 4
Proof of Lemma 4.1. We use backward induction on t to prove that ∇FtT (x) = Gt (x) for all x ∈ Rd and t ∈ [T ]. As for the base case corresponding to t = T , fix any x ∈ Rd , and observe that ∇FTT (x) = ∇ET −1 [fT (ξT , x)] = ET −1 [∇x fT (ξT , x)] = GT (x). Here, the first equality follows from the construction of FTT (x) in Definition C.1. Next, the gradient may be interchanged with the conditional expectation by the dominated convergence theorem, which applies thanks to Assumptions 3.2 and 3.3. The resulting expression coincides with the definition of GT (x). T (x) = G As for the induction step corresponding to any t ∈ [T − 1], assume that ∇Ft+1 t+1 (x) for d d all x ∈ R . Next, fix an arbitrary x ∈ R , and observe that T ∇FtT (x) =∇xt Et−1 [ft (ξt , Ft+1 (x))] T T =Et−1 ∇Ft+1 (x)∇xt ft (ξt , Ft+1 (x)) T =Et−1 Gt+1 (x)∇xt ft (ξt , Ft+1 (x)) = Gt (x).
Here, the first equality follows from the definition of FtT , while the second equality follows from the dominated convergence theorem and the chain rule. The third and the fourth equalities then exploit the induction hypothesis and the definition of Gt , respectively. The claim ultimately follows because F (x) = F1T (x). Proof of Lemma 4.5. Throughout the proof we fix x, and for simplicity we suppress the dependence on x t in all functions. We then find an upper bound on νt2 and a recursive formula for Et by backward induction b T and Assumption 4.4 imply that on t. As for the base case corresponding to t = T , the definition of G T b T ∥22T ≤ ν 2TT =: ET < ∞. νT2 = ET −1 ∥G As for the induction step corresponding to any t with t + 1 ∈ [T ], we introduce the auxiliary variables Ss :=
Ms X
1
2(2s−1 )(ρs +1)ℓ qs (ℓ)2s −1 ℓ=0
and
Ss′ :=
Ms X
1
22s (ρs−1 +1)ℓ qs (ℓ)2s (ρs−1 +1)−1 ℓ=0
∀s ∈ [T − 1], (38)
which allows us to express the induction hypothesis compactly as TY −1 2t+1 b t+1 ∥2t+1 ≤ Et+1 νt+1 = Et ∥G max Ss , Ss′ . 2 s=t+1
53
(39)
t
b t in (23) to obtain In order to bound νt2 , we use the law of total expectation and the definition of G t b t ∥22t νt2 = Et−1 ∥G
Mt X
1 λt ,e 1 λt ,o 2t 1 λt gb − gbt − gbt = qt (ℓ) Et−1 λt = ℓ 2 2 2 qt (λt )2t t ℓ=0 Mt X 1 ℓ,e 1 ℓ,o 2t 1 ℓ Et−1 gbt − gbt − gbt = . 2 2 2 qt (ℓ)2t −1
(40)
ℓ=0
In the second line we have eliminated the condition on the event λt = ℓ, which is allowed because λt is independent of all other random objects. Next, we define ∆ℓ := gbtℓ − 21 gbtℓ,e − 12 gbtℓ,o and bℓ H b ℓ,e H b ℓ,o H b t+1 − 1 ∇xt ft ξt , E b t+1 − 1 ∇xt ft ξt , E b t+1 Aℓt := ∇xt ft ξt , E t t t 2 2 to avoid clutter. By the construction of gbtℓ , gbtℓ,e and gbtℓ,o in Definition 4.2, we thus have bℓ G b ℓ,e G bℓ H b ℓ,e H b t+1 ∇xt ft ξt , E b t+1 − 1 E b t+1 ∇xt ft ξt , E b t+1 ∆ℓ = E t t t t 2 1 b ℓ,o b ℓ,o b b − Et Gt+1 ∇xt ft ξt , Et Ht+1 2 h i ℓ b bℓ G b b b b = E − E [ G ] ∇ f ξ , E H − ∇ f ξ , E [ H ] t+1 t t+1 xt t t t+1 xt t t t t+1 t t h i 1 b ℓ,e b b ℓ,e H b t+1 ] ∇xt ft ξt , E b t+1 − ∇xt ft ξt , Et [H b t+1 ] − Et Gt+1 − Et [G t 2 h i 1 b ℓ,o b b ℓ,o H b t+1 ] Aℓt , b t+1 ] ∇xt ft ξt , E b t+1 − ∇xt ft ξt , Et [H b t+1 ] + Et [G − Et Gt+1 − Et [G t 2 where the second equality is obtained by adding terms that sum to 0 P and using the antithetic averP4 sample 4 p p p−1 age identity (17). Next, we observe that Jensen’s inequality yields | i=1 ai | ≤ 4 i=1 |ai | for all a1 , . . . , a4 ∈ R and p ≥ 1. Hence, the (p = 2t )-th conditional moment of ∆ℓ obeys the estimate t Et−1 ∥∆ℓ ∥22 h h i 2t i t ℓ b bℓ G b b b b ≤ 42 −1 Et−1 E − E [ G ] ∇ f ξ , E H − ∇ f ξ , E [ H ] t+1 t t+1 x t t t+1 x t t t t+1 t t t t 2 h h t i 2 t i 2 −1 b ℓ,e G b ℓ,e H b t+1 − Et [G b t+1 ] ∇xt ft ξt , E b t+1 − ∇xt ft ξt , Et [H b t+1 ] + 4 2t Et−1 E t t 2 2 h h i 2t i t 2 −1 ℓ,o ℓ,o b G b H b t+1 − Et [G b t+1 ] ∇xt ft ξt , E b t+1 − ∇xt ft ξt , Et [H b t+1 ] + 4 2t Et−1 E t t 2 2 h i t t −1 ℓ 2 2 b +4 Et−1 Et [Gt+1 ] At 2 . To bound the first three terms above, we recall that 2ab ≤ a2 + b2 for all a, b ∈ R. Hence, we obtain t Et−1 ∥∆ℓ ∥22 h 2t+1 2t+1 i 2t −1 ℓ b bℓ G b b b b ≤ 4 2 Et−1 E −E [ G ] + ∇ f ξ , E H −∇ f ξ , E [ H ] t+1 t t+1 xt t t t+1 xt t t t t+1 t t 2 2 t+1 h t −1 2 2t+1 i 2 ℓ,e b b ℓ,e G b b t+1 −Et [G b t+1 ] b + 42t +1 Et−1 E + ∇ f ξ , E H −∇ f ξ , E [ H ] (41) xt t t t+1 xt t t t t+1 t t 2 2 2 h 2t+1 2t+1 i 2t −1 ℓ,o b b ℓ,o G b b t+1 −Et [G b t+1 ] b + 42t +1 Et−1 E f ξ , E [ H ] + ∇ f ξ , E H −∇ x t t t+1 x t t t t+1 t t t t 2 2 2 h i t t 2 −1 ℓ 2 b +4 Et−1 Et [Gt+1 ] At 2 . 54
The remainder of the proof proceeds in three steps, where we analyze separately the terms appearing in (41). b t+1 ] Aℓ , which appears on the last line of (41). In In Step (i), we bound the 2t -th conditional moment of Et [G t t+1 b t+1 . Step (ii), we bound the 2 -th conditional moments of the three gradient-difference terms involving H t+1 Finally, in Step (iii), we bound the 2 -th conditional moments of the three value-difference terms involvb t+1 . To complete the induction step, we then substitute all obtained bounds into (41) and (40). ing G Step (i)
b t+1 ] Aℓt , we first observe that To bound the 2t -th conditional moment of Et [G h i h i t t t b t+1 ] Aℓ 2 ≤ Et−1 Et [G b t+1 ] 2 Aℓ 2 , Et−1 Et [G t 2 t 2 2
(42)
where the inequality follows from the submultiplicativity of the Frobenius norm established in Lemma A.3. We emphasize that entry-wise matrix p-norms are not submultiplicative for p > 2. This is our reason for b t (x) using the Frobenius norm (corresponding to p = 2) in the definition of the p-th conditional moment of G t ℓ 2 in (25). Next, we bound ∥At ∥2 . Jensen’s inequality readily implies that t ∥Aℓt ∥22 =
dX t−1
2 Aℓt ei 2
i=1
dt−1
2t−1
2t−1 −1
≤ (dt−1 )
X i=1
t
(43)
∥Aℓt ei ∥22 ,
where ei stands for the i-th standard basis vector in Rdt−1 . By the definition of Aℓt , we can rewrite Aℓt ei as bℓ H b ℓ,e H b ℓ,o H b t+1 − 1 ∇xt ft,i ξt , E b t+1 b t+1 − 1 ∇xt ft,i ξt , E ∇xt ft,i ξt , E t t t 2 2 ℓ ℓ b 2 b b H b t+1 −Et [H b b t+1 ] b t+1 ] E =∇xt ft,i ξt , Et Ht+1 −∇xt ft,i ξt , Et [Ht+1 ] −∇xt ft,i ξt , Et [H t i ℓ,e 1h b H b ℓ,e H b t+1 −Et [H b t+1 ] b t+1 ] −∇xt ft,i ξt , E b t+1 −∇2x ft,i ξt , Et [H b t+1 ] E + ∇xt ft,i ξt , Et [H t t t 2 h i 1 b ℓ,o H b ℓ,o H b t+1 ] −∇xt ft,i ξt , E b t+1 −∇2 ft,i ξt , Et [H b t+1 ] E b t+1 −Et [H b t+1 ] , + ∇xt ft,i ξt , Et [H xt t t 2 where the equality follows from the antithetic sample average identity (17). The norm of Aℓt ei thus satisfies bℓ H b t+1 − ∇xt ft,i ξt , Et [H b t+1 ] ∥Aℓt ei ∥2 ≤ ∇xt ft,i ξt , E t ℓ b H b t+1 ] E b t+1 − Et [H b t+1 ] − ∇2x ft,i ξt , Et [H t
t
2
1 b ℓ,e H b t+1 ] − ∇xt ft,i ξt , E b t+1 + ∇xt ft,i ξt , Et [H t 2 ℓ,e b H b t+1 ] E b t+1 − Et [H b t+1 ] − ∇2xt ft,i ξt , Et [H t 2 1 ℓ,o b H b t+1 ] − ∇xt ft,i ξt , E b t+1 + ∇xt ft,i ξt , Et [H t 2 ℓ,o b H b t+1 ] E b t+1 − Et [H b t+1 ] − ∇2xt ft,i ξt , Et [H t 2
by virtue of the triangle inequality. Next, we bound the three norm terms by using the estimate (24), which is allowed because the Hessians are Hölder continuous thanks to Assumption 4.3. We therefore obtain ∥Aℓt ei ∥2 ≤
Rt Rt bℓ H b ℓ,e H b t+1 − Et [H b t+1 ] ρt +1 + b t+1 − Et [H b t+1 ] ρt +1 E E t t 2 2 ρt + 1 2(ρt + 1) Rt b ℓ,o H b t+1 − Et [H b t+1 ] ρt +1 . + E t 2 2(ρt + 1)
55
Two successive applications of Jensen’s inequality yield |a + b + c|p ≤ 2p−1 (|a|p + 2p−1 (|b|p + |c|p )) for all a, b, c ∈ R and p ≥ 1. By raising both sides of the previous inequality to the power p = 2t we thus find t
R 2t t h t t bℓ H b t+1 − Et [H b t+1 ] 2 (ρt +1) 22 −1 E t 2 ρt + 1 i t t 1 b ℓ,e b b ℓ,o H b t+1 ] 2 (ρt +1) + 1 E b t+1 − Et [H b t+1 ] 2 (ρt +1) + Et Ht+1 − Et [H t 2 2 2 2 h R 2t t t t−1 t bℓ H b t+1 − Et [H b t+1 ] 2t (ρt +1) 22 −1 (dt )2 (ρt +1)−1 E ≤ t 2 (ρt +1) ρt + 1 i t t 1 b ℓ,e b b ℓ,o H b t+1 ] 2t (ρt +1) + 1 E b t+1 − Et [H b t+1 ] 2t (ρt +1) , + Et Ht+1 − Et [H t 2 (ρt +1) 2 (ρt +1) 2 2
∥Aℓt ei ∥22 ≤
where the second inequality follows from Lemma A.2 with p = 2, q = 2t (ρt + 1), m = 1 and n = dt . t Substituting the above bound on ∥Aℓt ei ∥22 back into (43) then results in the estimate t
∥Aℓt ∥22 ≤
R 2t t t−1 t−1 t 22 −1 (dt−1 )2 (dt )2 (ρt +1)−1 Āℓt , ρt + 1
(44)
where we have used the shorthand t bℓ H b t+1 − Et [H b t+1 ] 2t (ρt +1) Āℓt = E t 2 (ρt +1) t t 1 b ℓ,e b b ℓ,o H b t+1 ] 2t (ρt +1) + 1 E b t+1 − Et [H b t+1 ] 2t (ρt +1) . + Et Ht+1 − Et [H t 2 (ρt +1) 2 (ρt +1) 2 2 Substituting (44) into (42) and applying the Cauchy-Schwarz and Jensen inequalities then yields i h t t b t+1 ] 2 Aℓ 2 Et−1 Et [G t 2 2 R 2t t i h t t−1 t−1 t b t+1 ] 2 Āℓ ≤ 22 −1 (dt−1 )2 (dt )2 (ρt +1)−1 Et−1 Et [G t 2 ρt + 1 1 R 2t t h h i t+1 2 i1 t−1 t−1 t ℓ 2 2 b t+1 2 ≤ E Ā . 22 −1 (dt−1 )2 (dt )2 (ρt +1)−1 Et−1 G t−1 t 2 ρt + 1
(45)
b t+1 and H b t+1 are constructed from the same sets of samples, which implies that Et [G b t+1 ] and Āℓ Note that G t are dependent. To further simplify the bound in (45), we use the Minkowski inequality to obtain h Et−1
i 12
2 Āℓt
i1 2t+1 (ρt +1) 2 bℓ H b b E − E [ H ] t+1 t t+1 2t (ρt +1) t h i1 h i1 t+1 t+1 1 b ℓ,e H b ℓ,o H b t+1 − Et [H b t+1 ] 2t (ρt +1) 2 + 1 Et−1 E b t+1 − Et [H b t+1 ] 2t (ρt +1) 2 + Et−1 E t t 2 (ρt +1) 2 (ρt +1) 2 2 1 h i t+1 1 bℓ H b t+1 − Et [H b t+1 ] 2t+1 (ρt +1) 2 ≤ (dt ) 2 Et−1 E t 2 (ρt +1) ≤ Et−1
h
h i1 1 h i1 1 2t+1 (ρt +1) 2 2t+1 (ρt +1) 2 ℓ,e b ℓ,o b b b b b + Et−1 Et Ht+1 − Et [Ht+1 ] 2t+1 (ρt +1) + Et−1 Et Ht+1 − Et [Ht+1 ] 2t+1 (ρt +1) . 2 2
56
Here, the second inequality follows from Lemma A.2 with p = 2t (ρt + 1), q = 2t+1 (ρt + 1), m = 1 and n = dt . By using Lemma A.4 with p = 2t+1 (ρt + 1) and n = 2ℓ or n = 2ℓ−1 , we further obtain 1
h
i1 ℓ 2 2
Et−1 Āt
≤ (dt )
B22t+1 (ρt +1)
1 2
2
t
1 2 2t+1 (ρt +1) t−1) (2 (ρt +1)(ℓ−1)
B
1
+ (dt ) 2
i1 t+1 b t+1 2t+1 (ρt +1) 2 H 2 (ρt +1)
h
Et−1 (2t−1 )(ρ +1)ℓ
2
i1 t+1 b t+1 2t+1 (ρt +1) 2 H 2 (ρt +1)
h Et−1 1
= 1+2
2t−1 (ρt +1)
(dt )
B22t+1 (ρt +1)
1 2
(2t−1 )(ρt +1)ℓ
2
h Et−1
i1 t+1 b t+1 2t+1 (ρt +1) 2 . H 2 (ρt +1)
Substituting the resulting expression into (45) and combining it with (42) finally yields h Et−1
t
b t+1 ]Aℓt 2 Et [G 2
i
≤
R 2t t h i1 t+1 2 t−1 t−1 t b t+1 2 22 −1 (dt−1 )2 (dt )2 (ρt +1)−1 Et−1 G 2 ρt + 1 1
× 1+2
2t−1 (ρt +1)
(dt )
1 2
B22t+1 (ρt +1) 2(2t−1 )(ρt +1)ℓ
h Et−1
i1 2t+1 (ρt +1) 2 b Ht+1 2t+1 (ρt +1) .
We have thus established an upper bound on the last term in (41), which depends on conditional moments b t+1 and H b t+1 and decays geometrically with ℓ. This concludes Step (i). of G Step (ii)
b t+1 , we note that To bound the gradient difference terms in (41) depending on H
h i t+1 2t+1 i t+1 bℓ H bℓ H b t+1 − ∇xt ft ξt , Et [H b t+1 ] b t+1 − Et [H b t+1 ] 2 ∇xt ft ξt , E ≤ St2 Et−1 E t t 2 2 h i h i t+1 t+1 B t+1 t t+1 t t+1 bℓ H b t+1 2t+1 , b t+1 − Et [H b t+1 ] 2t+1 ≤ St2 (dt )2 −1 2 t Et−1 H ≤ St2 (dt )2 −1 Et−1 E t 2 2 22 ℓ h
Et−1
where the first inequality exploits Assumption 3.3, and the second inequality follows from Lemma A.2 with p = 2, q = 2t+1 , m = 1 and n = dt . The third inequality follows from the law of total expectation and from Lemma A.4 with p = 2t+1 and n = 2ℓ . A similar reasoning implies that h i 2t+1 i 2t+1 2t+1 2t −1 B2t+1 b ℓ,e H b t+1 ] b t+1 − ∇xt ft ξt , Et [H b ∇xt ft ξt , E E H ≤ S (d ) t t−1 t+1 2t+1 , t t 2 22t (ℓ−1) i h h 2t+1 i 2t+1 2t+1 2t −1 B2t+1 b ℓ,o H b b t+1 − ∇xt ft ξt , Et [H b t+1 ] Et−1 ∇xt ft ξt , E E H ≤ S (d ) t−1 t+1 2t+1 . t t t 2 22t (ℓ−1) h
Et−1
As the empirical averages here are taken only over samples with even or odd indices, we set n = 2ℓ−1 when applying Lemma A.4. In summary, we have found upper bounds on the gradient difference terms in (41), b t+1 and decay geometrically with ℓ. This concludes Step (ii). which depend on conditional moments of H Step (iii)
b t+1 , we note that To bound the terms in (41) depending on G h Et−1
h 2t+1 i 2t+1 i 2t −1 ℓ b bℓ G b b t+1 − Et [G b t+1 ] b E − E [ G ] E ≤ (d d) E G t t−1 t+1 t t+1 t t 2 2t+1 h i t+1 B t t+1 2 b t+1 t+1 ≤ (dt d)2 −1 22t ℓ Et−1 G 2 2 h i 2t+1 2t −1 B2t+1 b ≤ (dt d) Et−1 Gt+1 2 , 22t ℓ 57
where the first inequality follows from Lemma A.2 with p = 2, q = 2t+1 , m = dt and n = d, while the second inequality follows from the law of total expectation and from Lemma A.4 with p = 2t+1 and n = 2ℓ . The third inequality follows from Jensen’s inequality (or from Lemma A.2 with p = 2 and q = 2t+1 ). A similar reasoning involving Lemma A.4 with p = 2t+1 and n = 2ℓ−1 implies that i i h h t+1 2t+1 2t −1 B2t+1 b ℓ,e G b t+1 − Et [G b t+1 ] 2 b ≤ (d d) G Et−1 E E , t t+1 2 t−1 t 2 22t (ℓ−1) i i h h t+1 2t+1 2t −1 B2t+1 b ℓ,o G b t+1 − Et [G b t+1 ] 2 b . Et−1 E G E ≤ (d d) t t+1 t−1 t 2 2 22t (ℓ−1) b t+1 and decay geometrically with ℓ, completing Step (iii). These bounds depend on conditional moments of G Having analyzed all terms in (41) in Steps (i), (ii) and (iii), we are now ready to combine the resulting estimates to a bound on the conditional moments of ∆ℓ . For notational convenience, we define the constants t+1 −3
Yt :=3 · 22
t+1 −3
Yt′ :=3 · 22
t
(dt d)2 −1 B2t+1 , t+1
t
St2 (dt )2 −1 B2t+1 , R 2t 1 1 t+1 t t−1 t−1 t−1 t Yt′′ :=22 +2 −3 (dt−1 )2 (dt )2 (ρt +1)− 2 1 + 22 (ρt +1) B22t+1 (ρt +1) , ρt + 1 which are independent of ℓ. Substituting the bounds from Steps (i), (ii) and (iii) into (41) then yields h i h i t+1 Yt Yt′ t 2t+1 b t+1 2 b Et−1 ∥∆ℓ ∥22 ≤ 2t ℓ Et−1 G + E H t−1 t+1 2t+1 2 2 22t ℓ (47) 1 h h i1 i Yt′′ 2t+1 (ρt +1) 2 2t+1 2 b b Et−1 Ht+1 2t+1 (ρt +1) . + (2t−1 )(ρ +1)ℓ Et−1 Gt+1 2 t 2 Next, we bound the three conditional expectations appearing in (47). By Lemma 3.19, under Assumption 3.3, and using the fact that 2s−1 (ρs + 1) ≤ 2s , the second conditional expectation satisfies TY −1 b t+1 ∥2t+1 Et−1 ∥H ≤ D t+1 t+1 2
X Ms
s=t+1
where T Dt+1 = µT2
ℓ=0
TY −1 s=t+1
1 sℓ 2 2 qs (ℓ)2s −1
3Ss 2
2s
≤ Dt+1
s d2s −1 B2s+1
TY −1 s=t+1
Ss ,
and Ss is defined in (38). By adapting Lemma 3.19 in the obvious way, one can further show that h Et
TY −1 i t+1 b t+1 2t+1 (ρt +1) ≤ D′ H Ss′ , t+1 2 (ρt +1) s=t+1
where 2T (ρ +1) ′ := µT T −1 Dt+1
TY −1 s=t+1
3Ss 2
2s (ρs−1 +1)
s d2s (ρs−1 +1)−1 B2s+1 (ρs−1 +1)
b t+1 in (47) and Ss′ is defined in (38). Moreover, by the induction hypothesis, the conditional moment of G satisfies the bound (39). Substituting all of these conditional moment bounds into (47) then yields T −1 T −1 Yt′ Dt+1 Y Yt Et+1 Y t ′ Et−1 ∥∆ℓ ∥22 ≤ max Ss S , S + s s 22t ℓ s=t+1 22t ℓ s=t+1 q ′ TY −1 TY −1 Yt′′ Et+1 Dt+1 1 1 Et + (2t−1 )(ρ +1)ℓ max Ss , Ss′ 2 Ss′ 2 ≤ (2t−1 )(ρ +1)ℓ max Ss , Ss′ , t t 2 2 s=t+1 s=t+1
58
where the second inequality follows from the trivial observations that 2t−1 (ρt + 1) ≤ 2t and that both Ss as well as Ss′ are dominated by max{Ss , Ss′ }. The inequality also exploits the (recursive) definition of Et as q ′ . Et := Yt Et+1 + Yt′ Dt+1 + Yt′′ Et+1 Dt+1 t
Replacing the conditional moment Et−1 [∥∆ℓ ∥22 ] in (40) with the above bound finally yields t Mt TY −1 TY −1 X Et−1 ∥∆ℓ ∥22 ′ 2t max S , S ≤ E max Ss , Ss′ . νt = ≤ E S s s t t t t −1 2 qt (ℓ) s=t s=t+1 ℓ=0
This completes the induction step. By their recursive definition and thanks to Assumption 4.4, the constants Et , t ∈ [T − 1], are finite whenever ET is finite. Hence, the claim follows. Proof of Lemma 4.6. We simplify notation by suppressing the dependence on x for all functions of x in this b t (x) obeys the estimate proof. In the following we will show that the conditional bias of the estimator G b t ] − Gt Et−1 [G 2 ≤
T −1 X
1 2 L[t:s−1] Ss Et−1 [νs+1 ]2
s=t
1
1 T −1 X L[s+1:k−1] Sk Et−1 [(µ2k+1 )2 ] 2 Et−1 [µ2s+1 ] 2 + 2Mk +1 2Ms /2 k=s+1
(48)
b = E[G b 1 ], the claim then follows immediately. We prove (48) by backward induction for all t ∈ [T ]. As E[G] b T constitutes an unbiased estimator for GT . on t. As for the base case corresponding to t = T , note that G Thus, the left hand side of (48) evaluates to 0, and the claim holds trivially (in fact, the right hand side of (48) vanishes, too). As for the induction step corresponding to any t with t + 1 ∈ [T ], assume that 1 T −1 T −1 2 ] 21 X X L[s+1:k−1] Sk Et [(µ2k+1 )2 ] 2 E [µ 1 t s+1 2 b t+1 ] − Gt+1 ≤ Et [G + . L[t+1:s−1] Ss Et [νs+1 ] 2 2 2Mk +1 2Ms /2 s=t+1
k=s+1
In addition, note that X Mt 1 ℓ,e 1 ℓ,o 1 λt ,e 1 λt ,o λt ℓ Et−1 gbt − gbt − gbt gbt − gbt − gbt = 2 2 2 2 ℓ=0 h i b Mt [ G b Mt [ H b t+1 ]∇xt ft ξt , E b t+1 ] , = Et−1 [b gtMt ] = Et−1 E t t
b t ] = Et−1 Et−1 [G
1 qt (λt )
(49)
b t in (23), and the second uses the law of total where the first equality follows from the definition of G expectation along with the independence of λt from all other random objects. The third equality holds because gbtℓ,e and gbtℓ,o share the distribution of gbtℓ−1 conditional on Ft for all ℓ ∈ [Mt ] and because gbt0,e = gbt0,o = 0. Finally, the fourth equality follows from Definition 4.2. Next, we introduce the shorthands h i T b Mt [G b Mt [H b t+1 ] ∇xt ft ξt , E b t+1 ] − ∇xt ft ξt , Ft+1 A1 := Et−1 E t t 2 h i T b A2 := Et−1 Et [Gt+1 ] − Gt+1 ∇xt ft ξt , Ft+1 2
to simplify notation. Then, we have h h i i b Mt [ G b Mt [H b t ] − Gt = Et−1 E b t+1 ]∇xt ft ξt , E b t+1 ] − Et−1 Gt+1 ∇xt ft ξt , F T Et−1 [G t+1 t t 2 2 h i Mt b Mt b Mt b T b b b ≤ Et−1 Et [Gt+1 ]∇xt ft ξt , Et [Ht+1 ] − Et [Gt+1 ]∇xt ft ξt , Ft+1 (50) 2 h i b Mt [ G b t+1 ]∇xt ft ξt , F T − Gt+1 ∇xt ft ξt , F T = A1 + A2 , + Et−1 E t+1 t+1 t 2
59
where the first equality exploits (49) and the definition of Gt in (22). The second equality follows from b Mt [·]] = Et [·]. Next, we separately bound the the law of total expectation and the elementary identity Et [E t terms A1 and A2 . As for A1 , we use the Cauchy-Schwarz and Jensen inequalities to obtain h h i i b Mt G b Mt [H b t+1 b t+1 ] − ∇xt ft ξt , F T A1 ≤ Et−1 E f ξ , E ∇ t t x t t+1 t t 2 2 i h b Mt [H b t+1 ] − F T b t+1 E ≤ St Et−1 G t+1 2 t 2 1 h i h i1 2 2 2 2 Mt b T b b ≤ St Et−1 Gt+1 2 Et−1 Et [Ht+1 ] − Ft+1 2 i 1 i1 i h h h 2 2 2 2 2 Mt b T b b b b = St Et−1 Gt+1 2 Et−1 Et [Ht+1 ] − Et [Ht+1 ] 2 + Et−1 Et [Ht+1 ] − Ft+1 2 . Here, the second inequality follows from the smoothness of the integrand ft , as stipulated in Assumption 3.3. b Mt [·]. This is permissible because the second norm Note that we have also removed the first sample average E t term does not depend on ξt+1 and because the underlying samples are independent and identically distributed conditional on the information available at time t − 1. The third inequality and the equality follow from the Cauchy-Schwarz inequality and from a standard bias-variance decomposition of the mean squared error. By Lemma C.6, the bias term admits the following upper bound. TX −1 h i L[t+1:s−1] Ss Et [µ2s+1 ] 2 2 T b Et−1 Et [Ht+1 ] − Ft+1 2 ≤ Et−1 2Ms +1 s=t+1
In addition, the variance term satisfies h i h i b Mt [ H b t+1 ] − Et [H b t+1 ] 2 = 2−Mt Et−1 H b t+1 − Et [H b t+1 ] 2 ≤ 2−Mt Et−1 [µ2 ], Et−1 E t+1 t 2 2 b Mt [H b t+1 ] equals the variance of H b t+1 where the equality holds because the variance of the sample average E t M t b divided by the sample size 2 . The inequality follows from the observation that the variance of Ht+1 is bounded below by its second moment and from the definition of µ2t+1 in (19). Substituting the bounds on these mean and the variance terms into the above bound on A1 then yields !1 TX −1 2 ] 2 h i1 2 ] 2 L S E [µ E [µ s t t−1 t+1 [t+1:s−1] s+1 b t+1 2 2 Et−1 + A1 ≤ St Et−1 G 2 2Ms +1 2 Mt s=t+1 ! TX −1 2 ] 2 12 2 ] 21 L S E [µ E [µ 1 s t t−1 [t+1:s−1] s+1 t+1 2 ≤ St Et−1 [νt+1 ] 2 Et−1 , + Mt /2 2Ms +1 2 s=t+1 1
2 in (25) and because (a2 +b2 ) 2 ≤ |a|+|b| where the second inequality holds because of the definition of νt+1 for all a, b ∈ R. The Minkowski and Jensen inequalities further imply that TX 1 L2 2 2 2 1 −1 T −1 X 2 L[t+1:s−1] Ss Et [µ2s+1 ] 2 2 [t+1:s−1] Ss Et [µs+1 ] Et−1 ≤ E t−1 2Ms +1 22Ms +2 s=t+1
s=t+1
≤
T −1 X s=t+1
1
L[t+1:s−1] Ss Et−1 [(µ2s+1 )2 ] 2 . 2Ms +1
Using this inequality to simplify that above bound on A1 yields 1 1 TX −1 L[t+1:s−1] Ss Et−1 [(µ2s+1 )2 ] 2 Et−1 [µ2t+1 ] 2 1 2 2 A1 ≤ St Et−1 [νt+1 ] + . 2Ms +1 2Mt /2 s=t+1
60
As for the term A2 , we use the Cauchy-Schwarz and Jensen inequalities to obtain h i T b t+1 ] − Gt+1 A2 ≤ Et−1 ∇xt ft ξt , Ft+1 E [ G t 2 2 i h b t+1 ] − Gt+1 ≤ Lt Et−1 Et [G 2 1 1 TX TX −1 −1 L[s+1:k−1] Sk Et [(µ2k+1 )2 ] 2 Et [µ2s+1 ] 2 1 2 L[t+1:s−1] Ss Et [νs+1 ] 2 + ≤ Lt Et−1 , 2Mk +1 2Ms /2 s=t+1
k=s+1
where the second inequality follows from Assumption 3.2, and the third inequality exploits the induction hypothesis. Jensen’s inequality together with the law of total expectations further yields A2 ≤ Lt
T −1 X
1 2 L[t+1:s−1] Ss Et−1 [νs+1 ]2
s=t+1
1
Et−1 [µ2s+1 ] 2 + 2Ms /2
TX −1 k=s+1
1
L[s+1:k−1] Sk Et−1 [(µ2k+1 )2 ] 2 2Mk +1
.
Combining the bounds on A1 and A2 with (50) and noting that L[t:t−1] = 1, we finally obtain b t ] − Gt Et−1 [G 2
1 L[t+1:s−1] Ss Et−1 [(µ2s+1 )2 ] 2 2Ms +1 s=t+1 1 1 TX T −1 −1 X L[s+1:k−1] Sk Et−1 [(µ2k+1 )2 ] 2 Et−1 [µ2s+1 ] 2 1 2 2 + Lt L[t+1:s−1] Ss Et−1 [νs+1 ] + 2Mk +1 2Ms /2 s=t+1 k=s+1 1 1 T −1 T −1 X X L[s+1:k−1] Sk Et−1 [(µ2k+1 )2 ] 2 Et−1 [µ2s+1 ] 2 1 2 . = + L[t:s−1] Ss Et−1 [νs+1 ] 2 Mk +1 Ms /2 2 2 s=t 1 2 ≤ St Et−1 [νt+1 ]2
1
Et−1 [µ2t+1 ] 2 + 2Mt /2
TX −1
k=s+1
This completes the induction step, and thus the claim follows. b Proof of Theorem 4.9. To bound the root √ mean squared error of G(x) by ϵ, it suffices to bound both the bias b and the standard deviation of G(x) by ϵ/ 2 for any x ∈ X . As for the bias, Lemma 4.6 implies that 1
T −1 2 (x)] 21 X L[t+1:s−1] Ss E[(µ2s+1 (x))2 ] 2 E[µ t+1 b E[G(x)] − G1 (x) 2 ≤ + 2Ms +1 2Mt /2 t=1 s=t+1 1 s+1 1 2t+1 2t+1 T −1 T −1 X X 1 2t+1 t+1 L[t+1:s−1] Ss E (µ2s+1 (x))2 2s+1 E µt+1 (x) ≤ L[t−1] St E νt+1 (x) 2 + , 2Ms +1 2Mt /2 T −1 X
1 2 L[t−1] St E[νt+1 (x)] 2
t=1
s=t+1
where the last expression follows from repeated application of Jensen’s inequality and from the definitions of the conditional moments µpt (x) and νtp (x) in (19) and (25), respectively. Lemma 4.5 further implies that TY −1 1 1 2t+1 t+1 2t+1 2 ≤ Et+1 max Ss , Ss′ E νt+1 (x)
1 2t+1
,
s=t+1
where the auxiliary variables Ss and Ss′ are defined as in (38) in the proof of Lemma 4.5. Thus, we have Ms X
Ms X
s
Zs2 −1 Ss = = . 2(2s−1 )(ρs +1)ℓ qs (ℓ)2s −1 2(2s−1 )(ρs +1)ℓ (1 − rs )(2s −1)ℓ ℓ=0 ℓ=0 1
61
Here, the second equality holds because qs (ℓ) = (1−rs )ℓ /Zs , and Zs := (1−(1−rs )Ms +1 )/rs is shorthand for the normalization constant of the truncated geometric distribution Geo(rs |Ms ). Note that Zs ≤ 1/rs . The first upper bound on rs in Assumption 4.8 thus ensures that Ss remains bounded as Ms tends to infinity. Consequently, there exists a constant S̄s ≥ Ss that is independent of Ms . Similarly, one can show that there exists a constant S̄s′ ≥ Ss′ that is independent of Ms . In summary, we may then conclude that TY −1 1 1 2t+1 t+1 2t+1 E νt+1 (x) 2 max S̄s , S̄s′ ≤ Et+1
1 2t+1
s=t+1
uniformly across all x ∈ X and Ms , s = t + 1, . . . , T − 1. As ρs ≤ 1, Lemma 3.19 further implies that t+1
µ2t+1 (x) ≤ Dt+1
TY −1 s=t+1
S̄s
uniformly across all x ∈ X and Ms , s = t + 1, . . . , T − 1. Combining the above estimates then yields b ∥E[G(x)] − G1 (x)∥2 ≤
T −1 X
1 t+1
2 L[t−1] St Et+1
t=1
s=t+1
×
−1 TY
1 2t+1
Dt+1
QT −1
1 2t+1
s=t+1 S̄s 2Mt /2
1 2t+1
max S̄s , S̄s′ +
T −1 X s=t+1
1 1 QT −1 2s 2s L[t+1:s−1] Ss Ds+1 s′ =s+1 S̄s′ . 2Ms +1
By construction, all constants in the resulting bound are independent of the truncation points. Since Ms ≥ 0, we can lower bound the denominator 2Ms +1 by 2Ms /2 and introduce auxiliary constants Wt′ ≥ 0 that are independent of the truncation points to obtain a simpler but relaxed bias bound of the form b ∥E[G(x)] − G1 (x)∥2 ≤
T −1 X
Wt′
2Mt /2 t=1
.
Note that Wt′ absorbs all dependencies of the bias on the Lipschitz, smoothness√ and Hölder constants, as well ′ as on the dimensions, moment √ bounds and rate parameters. Setting Wt = Wt 2(T − 1), it is now evident that the bias drops below ϵ/ 2 uniformly across all x ∈ X if Mt = ⌈2 log2 (Wt /ϵ)⌉ for all t ∈ [T − 1]. Next, we show that the variance of the MLMC gradientQ estimator is bounded by ϵ2 /2. From the first part −1 2 of the proof we know that ν1 (x) is bounded above by E1 Tt=1 max{S̄t , S̄t′ } uniformly across all x ∈ X . 2 Thus, both supx∈X ν1 (x) and n1 are finite. Consequently, we may conclude that b V(G(x)) ≤
ν12 (x) ϵ2 ≤ n1 2
uniformly across all x ∈ X , where the two inequalities follow from (26) and the definition of n1 . Finally, the expected sampling cost of the MLMC gradient estimator is given by b E[C(G(x))] = n1
Mt TY −1 X t=1 ℓ=0
ℓ
qt (ℓ)2 = n1
Mt TY −1 X t=1 ℓ=0
T −1
Y 1 − (2 − 2rt )Mt +1 2ℓ (1 − rt )ℓ = n1 , Zt Zt (2rt − 1) t=1
where the first equality can be proved as in Lemma 3.22. As 2−1 < rt < 1 by virtue of Assumption 4.8, the resulting expression remains bounded as Mt tends to infinity. In summary, the above reasoning implies that the root mean squared error of the MLMC gradient estimator is bounded by ϵ uniformly across all x ∈ X if n1 = O(ϵ−2 ) and Mt = O(log(ϵ−1 )) for all t ∈ [T − 1], b in which case the expected sampling cost satisfies E[C(G(x))] ≤ O(ϵ−2 ). Hence, the claim follows. 62
Proof of Theorem 4.11. Under the stated assumptions, Lemmas 4.5, 4.6, and 4.7, together with Theorem 4.9, are applicable. It is therefore easy to verify that the claim is a direct consequence of [41, Theorem C.1].
63