Augmenting Automatic Differentiation for a Single-Server Queue via the Leibniz Integral Rule

Augmenting Automatic Differentiation for a Single-Server Queue via the Leibniz Integral Rule

Introduction and Motivation

This work presents a new approach for constructing unbiased recursive estimators for higher-order derivatives of the expected waiting time in a G/G/1 single-server queue, specifically targeting performance measures where existing sample path methods—such as Infinitesimal Perturbation Analysis (IPA)—fail to provide unbiased estimators beyond first-order derivatives. The approach leverages the Leibniz integral rule in conjunction with the Lindley recursion, offering a systematic augmentation of automatic differentiation (AD) to stochastic queueing systems.

Estimating unbiased derivatives for expectations in queueing systems is central to sensitivity analysis, simulation optimization, and for underpinning gradient-based algorithms in reinforcement learning and simulation-based control. While first-order derivative estimation is well-developed via IPA, Likelihood Ratio Methods (LRM), and related approaches, the construction of unbiased higher-order estimators—especially from a single sample path—remains challenging and under-studied, particularly for discontinuous sample functionals prevalent in queueing phenomena.

Setting, Model, and Limitations of Classical Estimators

The paper formulates the estimation problem as follows: for a FCFS G/G/1 queue with service time distribution parameter $\theta$ (affecting only $S_1$ for notational and analytic simplicity), and assuming i.i.d. continuous interarrival and service times, the goal is to estimate the $n$th derivative of $E[W_i]$ with respect to $\theta$ via unbiased sample estimates computed recursively along a realization.

The Lindley equation $W_{i+1} = (W_i + S_i - A_i)^+$ governs the queue evolution. Standard sample path differentiation yields recursive IPA estimator formulas for higher derivatives, but these are generally biased beyond first order, due to non-differentiabilities (kinks) at transition points and the inadequacy of the naïve product rule in the presence of indicator and max operators on stochastic arguments.

This bias is clearly exemplified in cases such as location and scale parameter perturbations: beyond first order, the IPA estimator returns zero, while analytic calculation indicates nonzero higher derivatives (e.g., for $A_1 \sim U(0,1)$ and $S_1 = \theta$, $d^2 E[W_2] / d\theta^2 = 1$ for $\theta \in (0,1)$, but IPA yields zero).

Augmented Estimation via the Leibniz Integral Rule

To overcome the systematic IPA bias for higher-order derivatives, the study applies the Leibniz rule for differentiating parameterized integrals. The key technical step is to move the expectation (specifically, integration over the interarrival time $S_1$0) inside the differentiation, and then apply Leibniz's rule to account for the effect of variable integration limits and parameter dependence.

For second derivatives, the estimator for $S_1$1 is

$S_1$2

where $S_1$3 is the PDF of $S_1$4. The first term is the standard IPA (or AD-generated) recursive estimator, and the second term represents the necessary augmentation to eliminate the bias. This formulation is shown to be unbiased under regularity assumptions (dominance, continuity of $S_1$5), and is extended recursively up to third derivatives:

$S_1$6

The proofs crucially rely on repeated application of Leibniz's integral and product rules, enabled by the independence structure and continuity assumptions on arrival and service distributions. The resulting estimators can be implemented recursively along a sample path, with parameter derivatives of the service time distribution evaluated analytically or by AD.

The approach generalizes to higher orders. For deterministic service times, the estimators recover the analytic derivatives of the mean waiting time with respect to $S_1$7 and for M/G/1 and other standard models, they yield non-trivial (often closed-form) recursive estimators for arbitrary order derivatives.

Implications and Theoretical Context

The methodology bridges a fundamental gap between simulation-based gradient estimation (IPA, LRM, conditional Monte Carlo, measure-valued differentiation) and symbolic/automatic differentiation for discontinuous functionals. By exposing and correcting the bias in higher-order IPA for queueing performance measures, it provides actionable unbiased estimators for sensitivity and optimization routines upon which high-dimensional simulation-based stochastic optimization/learning algorithms depend.

The independence and continuity conditions imposed on interarrival and service time distributions are sufficient for technical requirements (dominated convergence, differentiability of the density at the integration limit), but the methodology hints at further generalization to broader classes of queueing models, discontinuous sample functions, and multidimensional parameterizations.

Numerically, the recursive formulas are modular and add little complexity relative to IPA, aside from density evaluations. The authors mention preliminary results showing empirically competitive performance for second derivative estimators, foreshadowing future extensive comparisons to LRM, weak derivative, smoothed PA, and augmented IPA approaches.

Practical and Theoretical Extensions

The formalism opens several directions, including:

• Extension to parameter vectors and arbitrary parametric dependence across all service times (not just $S_1$8), which is essential for gradient-based learning in reinforcement learning and simulation optimization.
• Application to other performance metrics, e.g., total system time, mean queue length, or statistics over busy periods.
• Incorporation of arrival process parameterizations, e.g., sensitivity to arrival rate or renewal structure.
• Generalization to more complex networks of queues, non-renewal processes, and discontinuous functionals of the sample path.

On the theoretical side, the explicit link between AD, IPA, and measure-valued differentiation for non-smooth functionals can inform future research on unbiased gradient estimation in stochastic networks, rare event simulation, and policy gradient algorithms for partially observed Markov decision processes.

Conclusion

The paper establishes a principled, recursive methodology for unbiased higher-order derivative estimation of expected waiting time in a G/G/1 queue, augmenting traditional AD and IPA estimators using the Leibniz integral rule. By resolving inherent bias in classical IPA for discontinuous sample functionals, the estimators extend the arsenal of tools available for gradient-based simulation optimization and learning. The framework sets a path for broader extension to complex queueing systems, higher-dimensional parameterizations, and new gradient-based algorithms in stochastic dynamic systems.