Deep Learning Method for Stationary Distribution of Reflected Brownian Motion

# Deep Learning Method for Stationary Distribution of Reflected Brownian Motion
Source: [https://arxiv.org/html/2607.08091](https://arxiv.org/html/2607.08091)
Zhanhao Zhang Operations Research and Information Engineering, Cornell University, zz564@cornell\.edu

###### Abstract

The stationary distribution of reflected Brownian motion \(RBM\) plays an important role in the analysis of high\-dimensional stochastic systems, yet closed\-form solutions are known only for a few special cases\. Computing important performance metrics, such as tail probabilities, is even more intractable, despite their practical relevance\. In this paper, we develop a deep learning approach that accurately and efficiently learns the Laplace transform of high\-dimensional RBMs based on the basic adjoint relationship \(BAR\)\. Our framework combines a careful design of the loss function, training data sampling procedure, and neural network architecture\. We evaluate the proposed method on RBM instances with known ground\-truth tail probabilities and demonstrate near\-perfect prediction in high\-dimensional settings, highlighting its potential as a general tool for analyzing stochastic systems beyond analytically tractable regimes\. Our code can be found at[https://github\.com/zhangz73/NN4MGF](https://github.com/zhangz73/NN4MGF)\.

## 1Introduction

Reflected Brownian motion \(RBM\) plays an important role in the analysis of multiclass queueing networks, where it often arises as a diffusion approximation under heavy traffic\. In such settings, the stationary distribution of an RBM provides useful approximations for the steady\-state behavior of the underlying queueing network\. Closed\-form expressions for stationary distributions are known only for a few special classes of RBMs\. This paper develops a scalable deep learning method for computing the Laplace transform of the stationary distribution of a high dimensional RBM\. The computed Laplace transforms can be used to estimate tail probabilities that serve as important performance metrics, such as tail latency for queueing networks\.

We consider add\-dimensional RBMZ=\{Z​\(t\),t≥0\}Z=\\\{Z\(t\),t\\geq 0\\\}associated with data\(Σ,μ,R\)\(\\Sigma,\\mu,R\), which satisfies the following equations:

Z​\(t\)=Z​\(0\)\+X​\(t\)\+R​Y​\(t\),t≥0,\\displaystyle Z\(t\)=Z\(0\)\+X\(t\)\+RY\(t\),\\qquad t\\geq 0,X=\{X​\(t\),t≥0\}​is ad\-dimensional Brownian motion with\\displaystyle X=\\\{X\(t\),t\\geq 0\\\}\\text\{ is a $d$\-dimensional Brownian motion with\}covariance matrixΣ\\Sigmaand driftμ\\mu,Y​\(0\)=0,Y​\(⋅\)​is non\-decreasing,\\displaystyle Y\(0\)=0,\\quad Y\(\\cdot\)\\text\{ is non\-decreasing\},∫0∞Zk​\(t\)​𝑑Yk​\(t\)=0,k=1,…,d\.\\displaystyle\\int\_\{0\}^\{\\infty\}Z\_\{k\}\(t\)dY\_\{k\}\(t\)=0,\\qquad k=1,\\dots,d\.Thed×dd\\times dmatrixRRis known as the reflection matrix\. We assume that the RBM is well defined and has a unique stationary distributionπ\\pi, which is satisfied, e\.g\., whenRRis anMM\-matrix andR−1​μ<0R^\{\-1\}\\mu<0\(Harrison\-Williams 1987\)\. Define the Laplace transform ofZZat steady state as

φ0​\(θ\)=𝔼π​\[e⟨−θ,Z​\(0\)⟩\],\\displaystyle\\varphi\_\{0\}\(\\theta\)=\\mathbb\{E\}\_\{\\pi\}\\left\[e^\{\\langle\-\\theta,Z\(0\)\\rangle\}\\right\],θ∈ℝ\+d\\displaystyle\\theta\\in\\mathbb\{R\}^\{d\}\_\{\+\}φk​\(θ\)=𝔼π​\[∫01e⟨−θ,Z​\(t\)⟩​𝑑Yk​\(t\)\],\\displaystyle\\varphi\_\{k\}\(\\theta\)=\\mathbb\{E\}\_\{\\pi\}\\left\[\\int\_\{0\}^\{1\}e^\{\\langle\-\\theta,Z\(t\)\\rangle\}dY\_\{k\}\(t\)\\right\],k=1,…,d,\\displaystyle k=1,\\dots,d,whereZ​\(0\)Z\(0\)follows the stationary distributionπ\\pi\. Lemma 1 of\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\]shows that the Laplace transformsφ\\varphiandφk\\varphi\_\{k\}are uniquely characterized by the Laplace version of basic adjoint relationship \(BAR\)

γ0​\(θ\)​φ0​\(θ\)=\\displaystyle\\gamma\_\{0\}\(\\theta\)\\varphi\_\{0\}\(\\theta\)=∑k=1dγk​\(θ\)​φk​\(θ\),\\displaystyle\\sum\_\{k=1\}^\{d\}\\gamma\_\{k\}\(\\theta\)\\varphi\_\{k\}\(\\theta\),\(1\)whereγ0​\(θ\)=−12​⟨θ,Σ​θ⟩\+⟨μ,θ⟩\\gamma\_\{0\}\(\\theta\)=\-\\frac\{1\}\{2\}\\langle\\theta,\\Sigma\\theta\\rangle\+\\langle\\mu,\\theta\\rangleandγk​\(θ\)=−⟨R\(k\),θ⟩\\gamma\_\{k\}\(\\theta\)=\-\\langle R^\{\(k\)\},\\theta\\rangle\. Here,R\(k\)R^\{\(k\)\}denotes thekk\-th column ofRR\. In\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\], the authors prove that the Laplace version BAR \([1](https://arxiv.org/html/2607.08091#S1.E1)\) is equivalent to the PDE version of the BAR that was first advanced in Harrison\-Williams \(1987\)\.

In contrast to\[[8](https://arxiv.org/html/2607.08091#bib.bib6)\], which develops an efficient approach for estimating steady\-state expectations of high\-dimensional RBMs, our work targets the Laplace transforms of their stationary distributions\. Since the Laplace transform can be numerically inverted to recover tail probabilities and, more broadly, the full stationary distribution\[[2](https://arxiv.org/html/2607.08091#bib.bib9)\], it provides a significantly richer characterization of system performance\. To the best of our knowledge, this is the first framework for estimating the Laplace transformφk​\(⋅\)\\varphi\_\{k\}\(\\cdot\)of high\-dimensional RBMs\. The main methodological contribution is to turn the BAR characterization into a scalable learning problem by combining a carefully designed loss function tailored to the BAR and the structural properties of Laplace transforms, a targeted training data sampling scheme, and a neural network architecture whose number of trainable parameters does not scale with the dimensiondd\. This yields an accurate and efficient computational tool for evaluating important performance metrics, such as tail probabilities and tail latencies, in large\-scale stochastic systems beyond analytically tractable regimes\.

##### Literature review\.

With the advancement of computational power, there has been growing interest in developing numerical methods for stochastic systems\. For instance,\[[17](https://arxiv.org/html/2607.08091#bib.bib26)\]propose a deep learning approach to compute convergence rates of Markov chains, and\[[18](https://arxiv.org/html/2607.08091#bib.bib24)\]extend this framework to estimate Lyapunov functions, solve Poisson equations, and approximate stationary distributions\. However, these methods are currently limited to low\-dimensional settings \(e\.g\., two dimensions\)\. Extending deep learning methods to high\-dimensional stochastic systems remains significantly more challenging and has motivated a growing line of research in stochastic control and diffusion approximations\.

Our work is closely related to this literature on applying deep learning to high\-dimensional stochastic systems\. In particular, the deep BSDE framework\[[14](https://arxiv.org/html/2607.08091#bib.bib12),[13](https://arxiv.org/html/2607.08091#bib.bib13)\]provides a powerful approach for solving high\-dimensional stochastic control problems by leveraging connections between parabolic PDEs and backward stochastic differential equations\. This framework has been successfully applied to a wide range of settings, including matching\[[5](https://arxiv.org/html/2607.08091#bib.bib20)\], scheduling\[[6](https://arxiv.org/html/2607.08091#bib.bib22)\], impulse control\[[7](https://arxiv.org/html/2607.08091#bib.bib19)\], drift control\[[4](https://arxiv.org/html/2607.08091#bib.bib21)\], and singular control\[[3](https://arxiv.org/html/2607.08091#bib.bib23)\]\. Our work is also motivated by stochastic systems, as the equation \([1](https://arxiv.org/html/2607.08091#S1.E1)\) characterizes the stationary behavior of reflected Brownian motions\. However, our setting is fundamentally different: unlike these approaches, which often exploit stochastic representations, such as the Feynman–Kac representation, to reformulate high\-dimensional PDEs as stochastic differential equations, \([1](https://arxiv.org/html/2607.08091#S1.E1)\) describes a stationary relationship and does not admit such a reformulation\.

Our work is also related to the extensive literature on applying deep learning to high\-dimensional partial differential equations \(PDEs\); see the survey by\[[22](https://arxiv.org/html/2607.08091#bib.bib11)\]\. In particular, prior work has studied PDEs arising from variational formulations using deep learning\[[9](https://arxiv.org/html/2607.08091#bib.bib16),[23](https://arxiv.org/html/2607.08091#bib.bib17),[19](https://arxiv.org/html/2607.08091#bib.bib18)\]\. Our approach constructs a loss function based on the squared residual of \([1](https://arxiv.org/html/2607.08091#S1.E1)\), which parallels the least\-squares formulations used in\[[9](https://arxiv.org/html/2607.08091#bib.bib16),[19](https://arxiv.org/html/2607.08091#bib.bib18)\]\. However, directly minimizing this residual is insufficient in our setting due to generalizability and numerical stability issues\. To address this, we design a structured loss function with additional regularization terms that enforce key properties of the Laplace transform\. In addition, we propose a tailored sampling scheme and a neural network architecture that scale effectively to high\-dimensional problems\.

## 2Deep learning approach

In this section, we approximateφk​\(⋅\)\\varphi\_\{k\}\(\\cdot\)fork=0,…,dk=0,\\dots,din \([1](https://arxiv.org/html/2607.08091#S1.E1)\) using feedforward neural networks in the complex domain, where \([1](https://arxiv.org/html/2607.08091#S1.E1)\) continues to hold by analytic extension\. Working in the complex domain is necessary because robust numerical inversion methods, such as the Talbot method\[[20](https://arxiv.org/html/2607.08091#bib.bib4)\]and related approaches\[[21](https://arxiv.org/html/2607.08091#bib.bib8),[1](https://arxiv.org/html/2607.08091#bib.bib7)\], evaluate the Laplace transform at complex arguments when computing tail probabilities\. Letθ=θRe\+i​θIm∈ℂd\\theta=\\theta^\{\\mathrm\{Re\}\}\+i\\theta^\{\\mathrm\{Im\}\}\\in\\mathbb\{C\}^\{d\}be add\-dimensional complex vector, withθRe,θIm∈ℝd\\theta^\{\\mathrm\{Re\}\},\\theta^\{\\mathrm\{Im\}\}\\in\\mathbb\{R\}^\{d\}\. We focus on accurate approximation over the bounded region

Θ:=\[θ¯Re,θ¯Re\]d×\[−θ¯Im,θ¯Im\]d\.\\displaystyle\\Theta:=\[\\underline\{\\theta\}^\{\\mathrm\{Re\}\},\\bar\{\\theta\}^\{\\mathrm\{Re\}\}\]^\{d\}\\times\[\-\\bar\{\\theta\}^\{\\mathrm\{Im\}\},\\bar\{\\theta\}^\{\\mathrm\{Im\}\}\]^\{d\}\.A naive approach is to use shallow feedforward neural networks to approximate the functionsφk​\(⋅\)\\varphi\_\{k\}\(\\cdot\), and then update the network parameters by minimizing the mean\-squared error between the left\-hand side and the right\-hand side of \([1](https://arxiv.org/html/2607.08091#S1.E1)\), where the training samples are drawn uniformly fromΘ\\Theta\. This straightforward framework, however, suffers from several fundamental difficulties\.

##### Poor generalization\.

The equation \([1](https://arxiv.org/html/2607.08091#S1.E1)\) characterizes the Laplace transform of the stationary distribution, which is analytic and monotone along the real axis\. However, the naive training framework only minimizes a finite\-sample residual of \([1](https://arxiv.org/html/2607.08091#S1.E1)\) and does not enforce either analyticity or monotonicity\. As a result, the neural network may fit the sampled training points well while violating these structural properties and deviating significantly at unseen inputs\.

##### Numerical stability\.

Because the Laplace transformsφk​\(θ\)\\varphi\_\{k\}\(\\theta\)vary exponentially withθ\\theta, their values can differ by many orders of magnitude within\[θ¯,θ¯\]d\[\\underline\{\\theta\},\\bar\{\\theta\}\]^\{d\}, leading to numerical precision issues whenθ\\thetais either small or large\.

##### Imbalanced training data\.

Asddgrows, samples drawn uniformly from high\-dimensional boxes for the real and imaginary parts ofθ\\thetararely fall in corner regions, where many coordinates are simultaneously close to their extreme values\. Instead, a typical sample contains a mixture of small and large coordinate values across dimensions\. Consequently, the naive sampling scheme under\-represents regimes in which the Laplace transform exhibits large magnitudes or strong oscillatory behavior, making it difficult for the training algorithm to learn these regions accurately\.

##### Poor scalability\.

A standard feedforward network takes the full vectorθ∈ℂd\\theta\\in\\mathbb\{C\}^\{d\}as input, so the size of its input layer grows linearly withdd\. As a result, the first layer contains parameters tied to individual input coordinates, which must be learned from samples that capture sufficient variability ofθ\\thetaacross alldddimensions\. Asddincreases, achieving such coverage requires substantially larger sample sizes, making the naive architecture less scalable for high\-dimensional RBMs\.

In the rest of this section, we will describe how we address these issues through a careful design of the loss function, training data sampling, and neural network architecture\.

### 2\.1Loss function

Instead of directly learning the Laplace transformsφk​\(⋅\)\\varphi\_\{k\}\(\\cdot\), we parameterize their logarithms using neural networks\. Specifically, the networks output functionsfk​\(⋅\)f\_\{k\}\(\\cdot\)such that

φk​\(θ\)=exp⁡\(fk​\(θ\)\),k=0,…,d\.\\displaystyle\\varphi\_\{k\}\(\\theta\)=\\exp\(f\_\{k\}\(\\theta\)\),\\qquad k=0,\\dots,d\.This log\-parameterization improves numerical stability and allows us to work with quantities of comparable scale when evaluating the BAR equation\.

For anyθ∈Θ\\theta\\in\\Theta, we define the loss function

ℒ​\(θ\):=\\displaystyle\\mathcal\{L\}\(\\theta\):=ℒBAR​\(θ\)\+λpair⋅ℒpair​\(θ\)\+λmono⋅ℒmono​\(θ\)\\displaystyle\\mathcal\{L\}\_\{\{\\rm BAR\}\}\(\\theta\)\+\\lambda\_\{\\rm pair\}\\cdot\\mathcal\{L\}\_\{\\rm pair\}\(\\theta\)\+\\lambda\_\{\\rm mono\}\\cdot\\mathcal\{L\}\_\{\\rm mono\}\(\\theta\)\+\\displaystyle\+λCR⋅ℒCR​\(θ\)\+λzero⋅ℒzero,\\displaystyle\\lambda\_\{\\rm CR\}\\cdot\\mathcal\{L\}\_\{\\rm CR\}\(\\theta\)\+\\lambda\_\{\\rm zero\}\\cdot\\mathcal\{L\}\_\{\\rm zero\},\(2\)whereλpair,λmono,λCR,λzero\>0\\lambda\_\{\\rm pair\},\\lambda\_\{\\rm mono\},\\lambda\_\{\\rm CR\},\\lambda\_\{\\rm zero\}\>0are hyperparameters that balance the relative importance of the different loss components\. The first term enforces the BAR equation, while the remaining terms incorporate structural properties that the Laplace transforms are known to satisfy\.

#### 2\.1\.1Normalized BAR error\.

The primary objective is to enforce the BAR equation \([1](https://arxiv.org/html/2607.08091#S1.E1)\)\. However, directly comparing the two sides of the equation can lead to numerical instability because the Laplace transforms may vary exponentially withθ\\theta\. To mitigate this issue, we introduce a normalization factor\.

For any inputθ∈Θ\\theta\\in\\Theta, we first compute the normalization factorν​\(θ\)\\nu\(\\theta\)as:

ν​\(θ\):=\\displaystyle\\nu\(\\theta\):=maxk=0,…,d⁡\{log⁡\|γk​\(θ\)\|\+fkR​e​\(θ\)\}\.\\displaystyle\\max\_\{k=0,\\dots,d\}\\left\\\{\\log\|\\gamma\_\{k\}\(\\theta\)\|\+f^\{Re\}\_\{k\}\(\\theta\)\\right\\\}\.This normalization rescales the two sides of the BAR equation to a comparable magnitude before evaluating their difference\. We then compute

κl​\(θ\)=\\displaystyle\\kappa\_\{l\}\(\\theta\)=γ0​\(θ\)⋅exp⁡\(f0​\(θ\)−ν​\(θ\)\),\\displaystyle\\gamma\_\{0\}\(\\theta\)\\cdot\\exp\\Big\(f\_\{0\}\(\\theta\)\-\\nu\(\\theta\)\\Big\),κr​\(θ\)=\\displaystyle\\kappa\_\{r\}\(\\theta\)=∑k=1dγk​\(θ\)⋅exp⁡\(fk​\(θ\)−ν​\(θ\)\)\.\\displaystyle\\sum\_\{k=1\}^\{d\}\\gamma\_\{k\}\(\\theta\)\\cdot\\exp\\Big\(f\_\{k\}\(\\theta\)\-\\nu\(\\theta\)\\Big\)\.Finally, the normalized BAR error is given as

ℒBAR​\(θ\):=\\displaystyle\\mathcal\{L\}\_\{\\rm BAR\}\(\\theta\):=\(\|κl​\(θ\)−κr​\(θ\)\|\|κl​\(θ\)\|\+\|κr​\(θ\)\|\+ϵ\)2,\\displaystyle\\left\(\\frac\{\|\\kappa\_\{l\}\(\\theta\)\-\\kappa\_\{r\}\(\\theta\)\|\}\{\|\\kappa\_\{l\}\(\\theta\)\|\+\|\\kappa\_\{r\}\(\\theta\)\|\+\\epsilon\}\\right\)^\{2\},\(3\)whereϵ\>0\\epsilon\>0is a very small constant to avoid division by zero\. This normalized error stabilizes training by comparing the relative discrepancy between the two sides of \([1](https://arxiv.org/html/2607.08091#S1.E1)\)\.

#### 2\.1\.2Pairwise consistency penalty\.

Next, we exploit structural constraints implied by the BAR equation by constructing points at which only one boundary term remains active\. Given inputθ∈Θ\\theta\\in\\Theta, we constructθ~\(1\),…,θ~\(d\)\\tilde\{\\theta\}^\{\(1\)\},\\dots,\\tilde\{\\theta\}^\{\(d\)\}that satisfiesR−k​θ~\(k\)=0R\_\{\-k\}\\tilde\{\\theta\}^\{\(k\)\}=0andθ~k\(k\)=θk\\tilde\{\\theta\}^\{\(k\)\}\_\{k\}=\\theta\_\{k\}, whereR−kR\_\{\-k\}denotes the matrix obtained fromRRby removing itskk\-th row\. By construction, these vectors lie on the boundary where only one reflection term remains active\. In particular, for anyk,k′=1,…,dk,k^\{\\prime\}=1,\\dots,dsuch thatk≠k′k\\neq k^\{\\prime\}, we haveγk′​\(θ~\(k\)\)=0\\gamma\_\{k^\{\\prime\}\}\(\\tilde\{\\theta\}^\{\(k\)\}\)=0\.

The original BAR equation relates the interior function only to the aggregate contribution of all boundary functions\. Consequently, each boundary function is constrained only indirectly through this aggregate residual, making it difficult to learn the individual boundary functions efficiently\. We therefore use the constructed points above to directly relate each boundary function to the interior function, providing an explicit constraint for each individual boundary function\. Substituting this relation into \([1](https://arxiv.org/html/2607.08091#S1.E1)\) and taking logarithms yields the following pairwise consistency penalty, which substantially accelerates convergence:

ℒpair​\(θ\):=∑k=1d\(log⁡γ0​\(θ~\(k\)\)\+f0​\(θ~\(k\)\)−log⁡γk​\(θ~\(k\)\)−fk​\(θ~\(k\)\)\)2\.\\displaystyle\\mathcal\{L\}\_\{\\rm pair\}\(\\theta\):=\\sum\_\{k=1\}^\{d\}\\Big\(\\log\\gamma\_\{0\}\(\\tilde\{\\theta\}^\{\(k\)\}\)\+f\_\{0\}\(\\tilde\{\\theta\}^\{\(k\)\}\)\-\\log\\gamma\_\{k\}\(\\tilde\{\\theta\}^\{\(k\)\}\)\-f\_\{k\}\(\\tilde\{\\theta\}^\{\(k\)\}\)\\Big\)^\{2\}\.\(4\)

#### 2\.1\.3Monotonicity penalty\.

The Laplace transforms are non\-increasing in the real part ofθ\\theta\. To enforce this property in the learned functions, we introduce a monotonicity penalty\. Given any inputθ∈Θ\\theta\\in\\Theta, we constructθ~:=R​e​\(θ\)\+0​i\\tilde\{\\theta\}:=Re\(\\theta\)\+0i, which lies on the real axis\.

ℒmono\(θ\):=∑k=0d\(∂fkRe​\(θ~\)∂Re​\(θ~\)\)\+\+λimg∑k=0d\(fkIm\(θ~\)\)2,\\displaystyle\\mathcal\{L\}\_\{\\rm mono\}\(\\theta\):=\\sum\_\{k=0\}^\{d\}\\left\(\\frac\{\\partial f^\{\\rm Re\}\_\{k\}\(\\tilde\{\\theta\}\)\}\{\\partial\\rm Re\(\\tilde\{\\theta\}\)\}\\right\)^\{\+\}\+\\lambda\_\{\\rm img\}\\sum\_\{k=0\}^\{d\}\(f^\{\\rm Im\}\_\{k\}\(\\tilde\{\\theta\}\)\)^\{2\},\(5\)whereλimg\>0\\lambda\_\{\\rm img\}\>0is a hyperparameter\. The first term penalizes violations of monotonicity, while the second enforces real\-valued outputs whenθ\\thetalies on the real axis\.

#### 2\.1\.4Cauchy\-Riemann penalty\.

Since the Laplace transforms are analytic functions, their real and imaginary parts must satisfy the Cauchy–Riemann equations\. We therefore include the following penalty to enforce approximate analyticity:

ℒCR​\(θ\):=∑k=0d\[\(∂fkRe​\(θ\)∂Re​\(θ\)−∂fkIm​\(θ\)∂Im​\(θ\)\)2\+\(∂fkRe​\(θ\)∂Im​\(θ\)\+∂fkIm​\(θ\)∂Re​\(θ\)\)2\]\.\\displaystyle\\mathcal\{L\}\_\{\\rm CR\}\(\\theta\):=\\sum\_\{k=0\}^\{d\}\\left\[\\left\(\\frac\{\\partial f^\{\\rm Re\}\_\{k\}\(\\theta\)\}\{\\partial\\rm Re\(\\theta\)\}\-\\frac\{\\partial f^\{\\rm Im\}\_\{k\}\(\\theta\)\}\{\\partial\\rm Im\(\\theta\)\}\\right\)^\{2\}\+\\left\(\\frac\{\\partial f^\{\\rm Re\}\_\{k\}\(\\theta\)\}\{\\partial\\rm Im\(\\theta\)\}\+\\frac\{\\partial f^\{\\rm Im\}\_\{k\}\(\\theta\)\}\{\\partial\\rm Re\(\\theta\)\}\\right\)^\{2\}\\right\]\.\(6\)

#### 2\.1\.5Zero\-anchoring penalty\.

Finally, we anchor the interior transform at the origin\. Sinceφ0​\(0\)=1\\varphi\_\{0\}\(0\)=1, we enforce this normalization through

ℒzero:=‖exp⁡\(f0​\(0\)\)−1‖22\.\\displaystyle\\mathcal\{L\}\_\{\\rm zero\}:=\\left\\\|\\exp\(f\_\{0\}\(0\)\)\-1\\right\\\|\_\{2\}^\{2\}\.\(7\)

### 2\.2Sampling strategy

At each gradient update for minimizing the loss function \([2](https://arxiv.org/html/2607.08091#S2.E2)\), we sample a batch ofNNdata points fromΘ\\Thetato compute the loss and its gradient with respect to the neural network parameters\. As discussed above, uniform sampling becomes increasingly ineffective in high dimensions because it rarely samples targeted corner regions\. These regions are important for capturing large\-magnitude Laplace values in the real domain and strong oscillatory behavior in the imaginary domain, both of which can significantly affect the accuracy of numerical inversion methods such as the Talbot method\. To improve coverage of these regions while maintaining support over the full domainΘ\\Theta, we adopt the following two\-stage sampling procedure\.

We sample each training data pointθ∈Θ\\theta\\in\\Thetain a batch as follows:

1. 1\.Step 1:We first sample scalar reference thresholdsθ∗Re∼Unif​\[θ¯Re,θ¯Re\]\\theta^\{\\mathrm\{Re\}\}\_\{\*\}\\sim\{\\rm Unif\}\[\\underline\{\\theta\}^\{\\mathrm\{Re\}\},\\bar\{\\theta\}^\{\\mathrm\{Re\}\}\]andθ∗Im∼Unif​\[0,θ¯Im\]\\theta^\{\\mathrm\{Im\}\}\_\{\*\}\\sim\{\\rm Unif\}\[0,\\bar\{\\theta\}^\{\\mathrm\{Im\}\}\]\. These thresholds determine the targeted corner regions for the real and imaginary parts\.
2. 2\.Step 2:Conditional on these thresholds, we sampleθRe∼Unif​\[θ¯Re,θ∗Re\]d\\theta^\{\\mathrm\{Re\}\}\\sim\{\\rm Unif\}\[\\underline\{\\theta\}^\{\\mathrm\{Re\}\},\\theta^\{\\mathrm\{Re\}\}\_\{\*\}\]^\{d\}andθIm∼Unif​\[\(\[−θ¯Im,−θ∗Im\]∪\[θ∗Im,θ¯Im\]\)d\]\\theta^\{\\mathrm\{Im\}\}\\sim\{\\rm Unif\}\\big\[\(\[\-\\bar\{\\theta\}^\{\\mathrm\{Im\}\},\-\\theta^\{\\mathrm\{Im\}\}\_\{\*\}\]\\cup\[\\theta^\{\\mathrm\{Im\}\}\_\{\*\},\\bar\{\\theta\}^\{\\mathrm\{Im\}\}\]\)^\{d\}\\big\], and constructθ=θRe\+i​θIm\\theta=\\theta^\{\\mathrm\{Re\}\}\+i\\theta^\{\\mathrm\{Im\}\}\. This conditional sampling scheme increases the probability of drawing points near the lower corner of the real domain and near the outer corners of the imaginary domain\.

We use this procedure to generate all training samples in each batch, with the batch sizeNNchosen independently of the dimensiondd\.

### 2\.3Neural network architecture

To address the scalability challenges discussed earlier, we design a neural network architecture whose size does not grow with the dimensiondd, which encodes each coordinate using a shared feature extractor and aggregates the resulting representations to produce the interior and boundary functions \(see Figure[1\(a\)](https://arxiv.org/html/2607.08091#S2.F1.sf1)\)\.

##### Shared coordinate encoder\.

Givenθ=\[θ1,…,θd\]∈Θ\\theta=\[\\theta\_\{1\},\\dots,\\theta\_\{d\}\]\\in\\Theta, we first apply a shared encoder to each coordinateθj\\theta\_\{j\}, which combines two components \(see Figure[1\(b\)](https://arxiv.org/html/2607.08091#S2.F1.sf2)\)\. The first component computes Fourier featuresη​\(θj\)\\eta\(\\theta\_\{j\}\)of the real and imaginary parts ofθj\\theta\_\{j\}, which improves the ability of the network to represent oscillatory patterns\. The second component maps the coordinate indexjjto an embedding vectorzjdimz\_\{j\}^\{\\rm dim\}, allowing the network to distinguish different coordinates while keeping the encoder weights shared across all dimensions\.

##### Interior network\.

For each coordinatej=1,…,dj=1,\\dots,d, the corresponding Fourier featuresη​\(θj\)\\eta\(\\theta\_\{j\}\)and the coordinate embedding vectorzjdimz\_\{j\}^\{\\rm dim\}are then passed to an interior network \(i\.e\. a feedforward neural network with two heads\)gintg\_\{\\rm int\}\. The network outputs a two\-dimensional vectorgint​\(η​\(θj\),zjdim\)=\(gintRe​\(η​\(θj\),zjdim\),gintIm​\(η​\(θj\),zjdim\)\)g\_\{\\rm int\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)=\\big\(g\_\{\\rm int\}^\{\\rm Re\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\),\\,g\_\{\\rm int\}^\{\\rm Im\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)\\big\), wheregintRe​\(η​\(θj\),zjdim\)g\_\{\\rm int\}^\{\\rm Re\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)\(resp\.gintIm​\(η​\(θj\),zjdim\)g\_\{\\rm int\}^\{\\rm Im\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)\) represents the real \(resp\. imaginary\) contribution of coordinatejjto the log of the interior Laplace transform\. These contributions are aggregated across coordinates to obtainf0Re​\(θ\)=∑j=1dgintRe​\(η​\(θj\),zjdim\)f^\{\\rm Re\}\_\{0\}\(\\theta\)=\\sum\_\{j=1\}^\{d\}g\_\{\\rm int\}^\{\\rm Re\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)andf0Im​\(θ\)=∑j=1dgintIm​\(η​\(θj\),zjdim\)f^\{\\rm Im\}\_\{0\}\(\\theta\)=\\sum\_\{j=1\}^\{d\}g\_\{\\rm int\}^\{\\rm Im\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\}\\big\)\. This additive structure allows the network to handle high\-dimensional inputs while maintaining a parameter count that does not scale withdd\.

##### Shared boundary network\.

To represent the boundary functions, we use another feedforward networkgbdaryg\_\{\\rm bdary\}with shared weights across all boundaries\. For each boundaryk=1,…,dk=1,\\dots,dand coordinatej=1,…,dj=1,\\dots,d, in addition to the Fourier featuresη​\(θj\)\\eta\(\\theta\_\{j\}\)and the coordinate embedding vectorzjdimz\_\{j\}^\{\\rm dim\}, the network also receives a boundary index embeddingzkbdaryz\_\{k\}^\{\\rm bdary\}\. The network outputs coordinate\-wise contributionsgbdary​\(η​\(θj\),zjdim,zkbdary\)g\_\{\\rm bdary\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\},z\_\{k\}^\{\\rm bdary\}\\big\), which are again aggregated to producefkRe​\(θ\)=∑j=1dgbdaryRe​\(η​\(θj\),zjdim,zkbdary\)f\_\{k\}^\{\\rm Re\}\(\\theta\)=\\sum\_\{j=1\}^\{d\}g\_\{\\rm bdary\}^\{\\rm Re\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\},z\_\{k\}^\{\\rm bdary\}\\big\)andfkIm​\(θ\)=∑j=1dgbdaryIm​\(η​\(θj\),zjdim,zkbdary\)f\_\{k\}^\{\\rm Im\}\(\\theta\)=\\sum\_\{j=1\}^\{d\}g\_\{\\rm bdary\}^\{\\rm Im\}\\big\(\\eta\(\\theta\_\{j\}\),z\_\{j\}^\{\\rm dim\},z\_\{k\}^\{\\rm bdary\}\\big\)fork=1,…,dk=1,\\ldots,d\.

Overall, this architecture leverages shared encoders and additive aggregation to achieve scalability in high dimensions while allowing the network to capture both interior and boundary behaviors of the Laplace transforms\.

![Refer to caption](https://arxiv.org/html/2607.08091v1/Plots/nn_pipeline.png)\(a\)Full neural network pipeline
![Refer to caption](https://arxiv.org/html/2607.08091v1/Plots/nn_coordinate_encoder.png)\(b\)Shared coordinate encoder

Figure 1:Neural network architecture

## 3Numerical Experiments

We evaluate the performance of our neural network by estimating tail probabilities of the formℙ​\(∑jZj\>t\)\\mathbb\{P\}\(\\sum\_\{j\}Z\_\{j\}\>t\)\. The predicted probabilities are obtained by applying numerical inverse Laplace transforms using the Talbot method\[[20](https://arxiv.org/html/2607.08091#bib.bib4)\]to the Laplace transform learned by the neural network\. To assess the accuracy of this approach, we consider two RBM examples for which ground truth results of tail probabilities can be obtained\.

The first example is a22\-dimensional RBM from\[[11](https://arxiv.org/html/2607.08091#bib.bib2)\]and\[[15](https://arxiv.org/html/2607.08091#bib.bib3)\], which admits a closed\-form expression for the stationary density but does not provide an explicit formula for the Laplace transform\. In this case, we compute the ground\-truth probabilitiesℙ​\(∑jZj\>t\)\\mathbb\{P\}\(\\sum\_\{j\}Z\_\{j\}\>t\)by numerically integrating the density function\. The second and third examples are2020\-dimensional and3030\-dimensional RBMs from\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\], whose Laplace transform admits a product\-form expression\. The corresponding ground\-truth probabilities are computed by applying the same Talbot inversion method to the exact Laplace transform\. Both numerical integration and Talbot inversion of Laplace transform can be achieved using thempmathpackage in Python\[[16](https://arxiv.org/html/2607.08091#bib.bib5)\]\.

![Refer to caption](https://arxiv.org/html/2607.08091v1/Plots/tail_prob_d=2_harrison.png)\(a\)22\-d RBM without closed\-form Laplace transform
![Refer to caption](https://arxiv.org/html/2607.08091v1/Plots/tail_prob_d=20.png)\(b\)2020\-d RBM with product form Laplace transform
![Refer to caption](https://arxiv.org/html/2607.08091v1/Plots/tail_prob_d=30.png)\(c\)3030\-d RBM with product form Laplace transform

Figure 2:Tail probabilities \(up to 1%\) of the stationary sum∑jZj\\sum\_\{j\}Z\_\{j\}: neural network vs\. ground truth\.The predicted tail probabilities match the ground truth almost perfectly in both examples \(see Figure[2](https://arxiv.org/html/2607.08091#S3.F2)\)\. The two\-dimensional example shows that our neural network can capture Laplace transforms with complex structures, where the behavior differs significantly across dimensions\. The2020\-dimensional and3030\-dimensional examples highlight the scalability of the architecture: despite the high dimensionality, the neural network still produces highly accurate probability estimates\. These results indicate that our neural network is both expressive enough to represent complex Laplace transforms and scalable to high\-dimensional RBMs\.

## 4Conclusion

In this paper, we propose a scalable deep learning framework for learning the Laplace transform of high\-dimensional RBMs\. Our approach integrates a carefully designed loss function, training data sampling procedure, and neural network architecture\. Numerical experiments demonstrate that, by numerically inverting the learned Laplace transform, our method can accurately compute tail probabilities of the stationary distribution in high\-dimensional settings\. These results highlight the potential of our framework as a general tool for estimating performance metrics in large\-scale stochastic systems\.

Despite these promising results, several limitations remain\. First, each gradient update currently requires16,38416\{,\}384training samples, leading to substantial GPU memory usage when processed in a single batch, or increased runtime when split across multiple batches, as the dimension of the RBM grows\. Developing more efficient training strategies to scale our approach to RBMs with hundreds or thousands of dimensions is an important direction for future work\. Second, beyond RBMs, it would be of interest to extend our framework to broader classes of stochastic systems and estimate practically relevant performance metrics, such as tail latencies and throughput\.

## Appendix AImplementation Details

##### Construction of Fourier features\.

To enhance the representation of complex inputs, we employ Fourier feature embeddings applied separately to the real and imaginary parts\. Consider an inputθ:=θRe\+i​θIm∈ℂd\\theta:=\\theta^\{\\rm Re\}\+i\\theta^\{\\rm Im\}\\in\\mathbb\{C\}^\{d\}\. For each coordinatej=1,…,dj=1,\\dots,d, we first normalize the inputs so that both the real and imaginary parts lie in\[−1,1\]\[\-1,1\]\. Specifically, we define the normalized real and imaginary componentsθ~jRe\\tilde\{\\theta\}^\{\\rm Re\}\_\{j\}andθ~jIm\\tilde\{\\theta\}^\{\\rm Im\}\_\{j\}as

θ~jRe:=\\displaystyle\\tilde\{\\theta\}^\{\\rm Re\}\_\{j\}:=2​θjRe−θ¯jReθ¯jRe−θ¯jRe−1,\\displaystyle 2\\frac\{\\theta^\{\\rm Re\}\_\{j\}\-\\underline\{\\theta\}^\{\\rm Re\}\_\{j\}\}\{\\bar\{\\theta\}^\{\\rm Re\}\_\{j\}\-\\underline\{\\theta\}^\{\\rm Re\}\_\{j\}\}\-1,θ~jIm:=\\displaystyle\\tilde\{\\theta\}^\{\\rm Im\}\_\{j\}:=θjIm\+θ¯jImθ¯jIm−1\.\\displaystyle\\frac\{\\theta^\{\\rm Im\}\_\{j\}\+\\bar\{\\theta\}^\{\\rm Im\}\_\{j\}\}\{\\bar\{\\theta\}^\{\\rm Im\}\_\{j\}\}\-1\.We then construct Fourier features based on these normalized inputs\. LetHHdenote the total number of Fourier features for each of the real and imaginary parts, whereHHis assumed to be divisible by22\. We construct two sets of frequency grids\{ξhRe\}h=1H/2\\\{\\xi\_\{h\}^\{\\rm Re\}\\\}\_\{h=1\}^\{H/2\}and\{ξhIm\}h=1H/2\\\{\\xi\_\{h\}^\{\\rm Im\}\\\}\_\{h=1\}^\{H/2\}as fixed, non\-trainable hyperparameters, using logarithmically spaced values over\[δminRe,δmaxRe\]\[\\delta\_\{\\min\}^\{\\rm Re\},\\delta\_\{\\max\}^\{\\rm Re\}\]and\[δminIm,δmaxIm\]\[\\delta\_\{\\min\}^\{\\rm Im\},\\delta\_\{\\max\}^\{\\rm Im\}\], respectively\. The resulting2​H2HFourier features are given by

\{sin⁡\(2​π​ξhRe​θjRe\),cos⁡\(2​π​ξhRe​θjRe\),sin⁡\(2​π​ξhIm​θjIm\),cos⁡\(2​π​ξhIm​θjIm\)\}h=1,…,H/2\.\\displaystyle\\Big\\\{\\sin\\Big\(2\\pi\\xi^\{\\rm Re\}\_\{h\}\\theta\_\{j\}^\{\\rm Re\}\\Big\),\\cos\\Big\(2\\pi\\xi^\{\\rm Re\}\_\{h\}\\theta\_\{j\}^\{\\rm Re\}\\Big\),\\sin\\Big\(2\\pi\\xi^\{\\rm Im\}\_\{h\}\\theta\_\{j\}^\{\\rm Im\}\\Big\),\\cos\\Big\(2\\pi\\xi^\{\\rm Im\}\_\{h\}\\theta\_\{j\}^\{\\rm Im\}\\Big\)\\Big\\\}\_\{h=1,\\dots,H/2\}\.

##### Hyperparameters for neural networks\.

In all experiments, we setH=64H=64, withδminRe=δminIm=−4\\delta\_\{\\min\}^\{\\rm Re\}=\\delta\_\{\\min\}^\{\\rm Im\}=\-4,δmaxRe=1\\delta\_\{\\max\}^\{\\rm Re\}=1, andδmaxIm=2\\delta\_\{\\max\}^\{\\rm Im\}=2\. The dimension index embeddingzjdimz\_\{j\}^\{\\rm dim\}and the boundary index embeddingzkbdaryz\_\{k\}^\{\\rm bdary\}are both chosen to have dimension6464\. For the interior network and the shared boundary network, we use feedforward architectures with two hidden layers, each consisting of128128neurons andSiLUactivation\.

##### Hyperparameters for training\.

For all experiments, we use theAdamWoptimizer with aCosineAnnealinglearning rate schedule, both implemented inPyTorch\. The learning rate is initialized at10−310^\{\-3\}and annealed to10−410^\{\-4\}\. We train the neural networks for100,000100\{,\}000epochs, where each epoch corresponds to a single gradient update\. For the2020\- and3030\-dimensional examples, we further fine\-tune the neural networks for an additional200,000200\{,\}000epochs\. Specifically, the learning rate is annealed from10−410^\{\-4\}to10−510^\{\-5\}over the first100,000100\{,\}000fine\-tuning epochs, and then from10−510^\{\-5\}to10−610^\{\-6\}over the remaining100,000100\{,\}000epochs\. At each epoch, we sample214=16,3842^\{14\}=16\{,\}384training data points to estimate the gradient of the loss function \([2](https://arxiv.org/html/2607.08091#S2.E2)\)\.

To reduce memory consumption, we evaluate the more expensive loss components \([4](https://arxiv.org/html/2607.08091#S2.E4)\)–\([6](https://arxiv.org/html/2607.08091#S2.E6)\) using random subsets of these16,38416\{,\}384training data points\. The pairwise consistency loss \([4](https://arxiv.org/html/2607.08091#S2.E4)\) is evaluated using only300300points\. This is because PyTorch stores the intermediate activations of the neural network for each input until backpropagation, and each inputθ\\thetais transformed intoddauxiliary inputsθ~\(1\),…,θ~\(d\)\\tilde\{\\theta\}^\{\(1\)\},\\dots,\\tilde\{\\theta\}^\{\(d\)\}, causing the memory requirement to grow linearly withdd\. For the monotonicity penalty \([5](https://arxiv.org/html/2607.08091#S2.E5)\) and Cauchy–Riemann penalty \([6](https://arxiv.org/html/2607.08091#S2.E6)\), we first use10241024points to evaluate the interior terms \(i\.e\.,k=0k=0\), and then use128128points to evaluate these penalties across allk=0,…,dk=0,\\dots,d\. These penalties are particularly memory intensive because they require first\-order derivatives of the neural network outputs, which require retaining the computational graph so that gradients can be backpropagated through these derivatives\. The penalty coefficients are set toλbdary=λmono=λCR=10\\lambda\_\{\\rm bdary\}=\\lambda\_\{\\rm mono\}=\\lambda\_\{\\rm CR\}=10andλzero=0\.1\\lambda\_\{\\rm zero\}=0\.1\.

##### Problem instances\.

We consider two classes of RBMs for all our numerical experiments\. The first is adopted from\[[11](https://arxiv.org/html/2607.08091#bib.bib2)\], where a closed\-form expression is available for the stationary distribution, but not for its Laplace transform or tail probabilities\. Specifically, we consider a22\-dimensional RBM with

Σ=\[1001\],μ=\[−10\],R=\[10−11\]\.\\displaystyle\\Sigma=\\begin\{bmatrix\}1&0\\\\ 0&1\\\\ \\end\{bmatrix\},\\qquad\\mu=\\begin\{bmatrix\}\-1\\\\ 0\\end\{bmatrix\},\\qquad R=\\begin\{bmatrix\}1&0\\\\ \-1&1\\end\{bmatrix\}\.The densityρ\\rhoof its stationary distribution is given by

ρ​\(θ\)=\\displaystyle\\rho\(\\theta\)=C​r1/2​e−\(r\+θ1\)​cos⁡\(ψ/2\),\\displaystyle Cr^\{1/2\}e^\{\-\(r\+\\theta\_\{1\}\)\}\\cos\(\\psi/2\),whereCC,rr, andψ\\psiare defined as

C=\\displaystyle C=1π​23/2,\\displaystyle\\frac\{1\}\{\\sqrt\{\\pi\}\}2^\{3/2\},r=\\displaystyle r=θ12\+θ22,\\displaystyle\\sqrt\{\\theta\_\{1\}^\{2\}\+\\theta\_\{2\}^\{2\}\},ψ=\\displaystyle\\psi=arccos⁡\(θ1/r\)\.\\displaystyle\\arccos\(\\theta\_\{1\}/r\)\.The tail probabilities are then obtained via numerical integration ofρ\\rhousingmpmath\.

The second class is adopted from\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\], which admits a product\-form Laplace transform due to theskew\-symmetryproperty\. An RBM is said to satisfy the skew\-symmetry condition if

2​Σ=R​diag​\(R\)−1​diag​\(Σ\)\+diag​\(Σ\)​diag​\(R\)−1​RT\.\\displaystyle 2\\Sigma=R\{\\rm diag\}\(R\)^\{\-1\}\{\\rm diag\}\(\\Sigma\)\+\{\\rm diag\}\(\\Sigma\)\{\\rm diag\}\(R\)^\{\-1\}R^\{T\}\.Following\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\], we construct an RBM for any dimensiond≥2d\\geq 2that satisfies the skew symmetry condition as follows:

Rj,j=1,\\displaystyle R\_\{j,j\}=1,∀j=1,…,d,\\displaystyle\\forall j=1,\\dots,d,Rj,j−1=−1,\\displaystyle R\_\{j,j\-1\}=\-1,∀j=2,…,d,\\displaystyle\\forall j=2,\\dots,d,Σj,j=cj\+cj\+1,\\displaystyle\\Sigma\_\{j,j\}=c\_\{j\}\+c\_\{j\+1\},∀j=1,…,d,\\displaystyle\\forall j=1,\\dots,d,Σj,j−1=−cj,\\displaystyle\\Sigma\_\{j,j\-1\}=\-c\_\{j\},∀j=2,…,d,\\displaystyle\\forall j=2,\\dots,d,Σj−1,j=−cj,\\displaystyle\\Sigma\_\{j\-1,j\}=\-c\_\{j\},∀j=2,…,d,\\displaystyle\\forall j=2,\\dots,d,μj=βj−βj\+1,\\displaystyle\\mu\_\{j\}=\\beta\_\{j\}\-\\beta\_\{j\+1\},∀j=1,…,d,\\displaystyle\\forall j=1,\\dots,d,where the vectorsc∈ℝd\+1c\\in\\mathbb\{R\}^\{d\+1\}andβ∈ℝd\+1\\beta\\in\\mathbb\{R\}^\{d\+1\}are defined by

cj=\\displaystyle c\_\{j\}=1,\\displaystyle 1,∀j=1,…,d\+1,\\displaystyle\\forall j=1,\\dots,d\+1,βj=\\displaystyle\\beta\_\{j\}=j,\\displaystyle j,∀j=1,…,d\+1\.\\displaystyle\\forall j=1,\\dots,d\+1\.According to\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\], the corresponding Laplace transforms are given by

φ0​\(θ\)=\\displaystyle\\varphi\_\{0\}\(\\theta\)=∏j=1dαjαj−θj,\\displaystyle\\prod\_\{j=1\}^\{d\}\\frac\{\\alpha\_\{j\}\}\{\\alpha\_\{j\}\-\\theta\_\{j\}\},∀θ<α,\\displaystyle\\forall\\theta<\\alpha,φk​\(θ\)=\\displaystyle\\varphi\_\{k\}\(\\theta\)=Σk,k2​Rk,k​αk​∏j≠kαjαj−θj,\\displaystyle\\frac\{\\Sigma\_\{k,k\}\}\{2R\_\{k,k\}\}\\alpha\_\{k\}\\prod\_\{j\\neq k\}\\frac\{\\alpha\_\{j\}\}\{\\alpha\_\{j\}\-\\theta\_\{j\}\},∀θ<α,\\displaystyle\\forall\\theta<\\alpha,whereα\\alphais defined as

α=−2​d​i​a​g​\(Σ\)−1​diag​\(R\)​R−1​μ\.\\displaystyle\\alpha=\-2\{\\rm diag\}\(\\Sigma\)^\{\-1\}\{\\rm diag\}\(R\)R^\{\-1\}\\mu\.The tail probabilities are then computed via numerical Laplace inversion using the Talbot method implemented inmpmath\.

## Appendix BMoments Estimation

### B\.1Numerical procedure for computing moments given Laplace transform

Given the Laplace transform, we compute the moments of the RBM stationary distribution by following the procedure in\[[10](https://arxiv.org/html/2607.08091#bib.bib10)\]\. To make this paper self\-contained, we summarize the key steps below\.

The method is based on numerical inversion of the Laplace transform in the complex plane\. For each moment ordernn, the inversion is carried out along a circular contour

𝒞n:=\{z∈ℂ:\|z\|=rn\},\\mathcal\{C\}\_\{n\}:=\\\{z\\in\\mathbb\{C\}:\|z\|=r\_\{n\}\\\},centered at the origin with radiusrnr\_\{n\}\. The contour integral representation of the moment is then approximated by evaluating the integrand at equally spaced points

zj=rn​ei​π​j/\(n​ℓ\),j=0,1,…,n​ℓ,z\_\{j\}=r\_\{n\}e^\{i\\pi j/\(n\\ell\)\},\\quad j=0,1,\\dots,n\\ell,which leads to a trapezoidal\-type discretization of the integral\. This yields the following inversion formula:

mn=n\!2​n​ℓ​rnn​ann​\[Wn​\(rn\)\+\(−1\)n​Wn​\(−rn\)\+2​∑j=1n​ℓ−1Re​\(Wn​\(zj\)​e−i​π​j\)\],\\displaystyle m\_\{n\}=\\frac\{n\!\}\{2n\\ell\\,r\_\{n\}^\{n\}\\,a\_\{n\}^\{n\}\}\\left\[W\_\{n\}\(r\_\{n\}\)\+\(\-1\)^\{n\}W\_\{n\}\(\-r\_\{n\}\)\+2\\sum\_\{j=1\}^\{n\\ell\-1\}\\rm Re\\\!\\left\(W\_\{n\}\\\!\(z\_\{j\}\)e^\{\-i\\pi j\}\\right\)\\right\],\(8\)wherern=10−ϵ/\(2​n​ℓ\)r\_\{n\}=10^\{\-\\,\\epsilon/\(2n\\ell\)\}for some accuracy parameterϵ\>0\\epsilon\>0and inversion parameterℓ∈\{1,2\}\\ell\\in\\\{1,2\\\}\. The parameterrnr\_\{n\}is chosen to control the discretization error arising from approximating the contour integral using a finite trapezoidal sum, while the parameterℓ\\ellhelps mitigate round\-off error by increasing the number of quadrature points along the contour\. The functionWn​\(z\):=φ​\(−an​z\)W\_\{n\}\(z\):=\\varphi\(\-a\_\{n\}z\)is a rescaled version of the Laplace transform\. The scaling factors\{an\}\\\{a\_\{n\}\\\}are chosen adaptively so that the magnitude ofWn​\(z\)W\_\{n\}\(z\)remains well\-conditioned on the contour𝒞n\\mathcal\{C\}\_\{n\}, preventing numerical overflow or underflow when computing high\-order moments\.

The overall procedure is summarized in Algorithm[1](https://arxiv.org/html/2607.08091#algorithm1)\. It first computes the first two moments, refines them using improved scaling, and then proceeds recursively to higher\-order moments using previously computed values to update the scaling factors\.

Input :Laplace transform

φ​\(θ\)\\varphi\(\\theta\), number of moments

M≥2M\\geq 2, accuracy parameter

ϵ\\epsilon, inversion parameter

ℓ∈\{1,2\}\\ell\\in\\\{1,2\\\}
// Step 1: Initial computation of

m1m\_\{1\}and

m2m\_\{2\}
Set

a1←1a\_\{1\}\\leftarrow 1and

a2←m1a\_\{2\}\\leftarrow m\_\{1\}
Compute

m1m\_\{1\}and

m2m\_\{2\}using \([8](https://arxiv.org/html/2607.08091#A2.E8)\) with

W1​\(z\)=φ​\(−a1​z\)W\_\{1\}\(z\)=\\varphi\(\-a\_\{1\}z\)and

W2​\(z\)=φ​\(−a2​z\)W\_\{2\}\(z\)=\\varphi\(\-a\_\{2\}z\)
// Step 2: Refinement of

m1m\_\{1\}and

m2m\_\{2\}
Set

a1←2​m1m2a\_\{1\}\\leftarrow\\frac\{2m\_\{1\}\}\{m\_\{2\}\}and

a2←2​m1m2a\_\{2\}\\leftarrow\\frac\{2m\_\{1\}\}\{m\_\{2\}\}
Recompute

m1,m2m\_\{1\},m\_\{2\}using \([8](https://arxiv.org/html/2607.08091#A2.E8)\) with

W1​\(z\)=φ​\(−a1​z\)W\_\{1\}\(z\)=\\varphi\(\-a\_\{1\}z\)and

W2​\(z\)=φ​\(−a2​z\)W\_\{2\}\(z\)=\\varphi\(\-a\_\{2\}z\)
// Step 3: Compute

mnm\_\{n\}for

3≤n≤M3\\leq n\\leq M
for*n=3,4,…,Mn=3,4,\\dots,M*do

Set

an←\(n−1\)​mn−2mn−1a\_\{n\}\\leftarrow\\frac\{\(n\-1\)m\_\{n\-2\}\}\{m\_\{n\-1\}\}
Compute

mnm\_\{n\}using \([8](https://arxiv.org/html/2607.08091#A2.E8)\) with

Wn​\(z\)=φ​\(−an​z\)W\_\{n\}\(z\)=\\varphi\(\-a\_\{n\}z\)
return

m1,m2,…,mMm\_\{1\},m\_\{2\},\\dots,m\_\{M\}

Algorithm 1Adaptive moment computation from Laplace transform
### B\.2Experiment results

Using the procedure described in Algorithm[1](https://arxiv.org/html/2607.08091#algorithm1), we estimate the moments of RBMs adopted from\[[12](https://arxiv.org/html/2607.08091#bib.bib1)\]with dimensions55,2020, and3030\. The results are reported in Tables[1](https://arxiv.org/html/2607.08091#A2.T1)–LABEL:tab:moments\-dim\-30\. Overall, the proposed framework provides accurate estimates for low\-order moments, particularly in lower\-dimensional settings\. For the55\-dimensional RBM, the estimated first\-, second\-, and third\-order moments all exhibit relatively small relative errors\. For the2020\- and3030\-dimensional RBMs, the first\-order moments also remain reasonably accurate across most coordinates\. As the dimension and moment order increase, however, the estimation errors become larger, especially for certain second\-order moments in the3030\-dimensional case\. This behavior is expected, as moment estimation relies on accurate local evaluations of the Laplace transform nearθ=0\\theta=0through contour\-based derivative approximation, which becomes increasingly sensitive to approximation errors in high\-dimensional settings\. In contrast, the Talbot inversion used for tail probability estimation tends to remain more stable in high dimensions because it depends on evaluating the Laplace transform over a broader region of the complex plane rather than accurately recovering local derivative information nearθ=0\\theta=0\. As a result, small local approximation errors are less significantly amplified compared with moment estimation\.

First MomentsDimensionjjTruePredAbs ErrRel Err01\.00×1001\.00\\times 10^\{0\}9\.98×10−19\.98\\times 10^\{\-1\}2\.06×10−32\.06\\times 10^\{\-3\}2\.06×10−32\.06\\times 10^\{\-3\}15\.00×10−15\.00\\times 10^\{\-1\}4\.99×10−14\.99\\times 10^\{\-1\}1\.19×10−31\.19\\times 10^\{\-3\}2\.38×10−32\.38\\times 10^\{\-3\}23\.33×10−13\.33\\times 10^\{\-1\}3\.34×10−13\.34\\times 10^\{\-1\}6\.46×10−46\.46\\times 10^\{\-4\}1\.94×10−31\.94\\times 10^\{\-3\}32\.50×10−12\.50\\times 10^\{\-1\}2\.51×10−12\.51\\times 10^\{\-1\}6\.83×10−46\.83\\times 10^\{\-4\}2\.73×10−32\.73\\times 10^\{\-3\}42\.00×10−12\.00\\times 10^\{\-1\}2\.00×10−12\.00\\times 10^\{\-1\}5\.72×10−55\.72\\times 10^\{\-5\}2\.86×10−42\.86\\times 10^\{\-4\}Second Moments02\.00×1002\.00\\times 10^\{0\}2\.00×1002\.00\\times 10^\{0\}3\.09×10−33\.09\\times 10^\{\-3\}1\.55×10−31\.55\\times 10^\{\-3\}15\.00×10−15\.00\\times 10^\{\-1\}4\.99×10−14\.99\\times 10^\{\-1\}8\.76×10−48\.76\\times 10^\{\-4\}1\.75×10−31\.75\\times 10^\{\-3\}22\.22×10−12\.22\\times 10^\{\-1\}2\.25×10−12\.25\\times 10^\{\-1\}2\.28×10−32\.28\\times 10^\{\-3\}1\.03×10−21\.03\\times 10^\{\-2\}31\.25×10−11\.25\\times 10^\{\-1\}1\.27×10−11\.27\\times 10^\{\-1\}2\.18×10−32\.18\\times 10^\{\-3\}1\.74×10−21\.74\\times 10^\{\-2\}48\.00×10−28\.00\\times 10^\{\-2\}7\.85×10−27\.85\\times 10^\{\-2\}1\.45×10−31\.45\\times 10^\{\-3\}1\.82×10−21\.82\\times 10^\{\-2\}Third Moments06\.00×1006\.00\\times 10^\{0\}5\.88×1005\.88\\times 10^\{0\}1\.16×10−11\.16\\times 10^\{\-1\}1\.93×10−21\.93\\times 10^\{\-2\}17\.50×10−17\.50\\times 10^\{\-1\}7\.27×10−17\.27\\times 10^\{\-1\}2\.31×10−22\.31\\times 10^\{\-2\}3\.08×10−23\.08\\times 10^\{\-2\}22\.22×10−12\.22\\times 10^\{\-1\}2\.27×10−12\.27\\times 10^\{\-1\}4\.71×10−34\.71\\times 10^\{\-3\}2\.12×10−22\.12\\times 10^\{\-2\}39\.38×10−29\.38\\times 10^\{\-2\}9\.60×10−29\.60\\times 10^\{\-2\}2\.22×10−32\.22\\times 10^\{\-3\}2\.36×10−22\.36\\times 10^\{\-2\}44\.80×10−24\.80\\times 10^\{\-2\}4\.32×10−24\.32\\times 10^\{\-2\}4\.84×10−34\.84\\times 10^\{\-3\}1\.01×10−11\.01\\times 10^\{\-1\}Table 1:Moment estimates ford=5d=5Table 2:Moment estimates ford=20d=20Table 3:Moment estimates ford=30d=30First MomentsDimensionjjTruePredAbs ErrRel Err01\.00×1001\.00\\times 10^\{0\}9\.90×10−19\.90\\times 10^\{\-1\}1\.04×10−21\.04\\times 10^\{\-2\}1\.04×10−21\.04\\times 10^\{\-2\}15\.00×10−15\.00\\times 10^\{\-1\}4\.75×10−14\.75\\times 10^\{\-1\}2\.48×10−22\.48\\times 10^\{\-2\}4\.97×10−24\.97\\times 10^\{\-2\}23\.33×10−13\.33\\times 10^\{\-1\}3\.19×10−13\.19\\times 10^\{\-1\}1\.47×10−21\.47\\times 10^\{\-2\}4\.42×10−24\.42\\times 10^\{\-2\}32\.50×10−12\.50\\times 10^\{\-1\}2\.37×10−12\.37\\times 10^\{\-1\}1\.32×10−21\.32\\times 10^\{\-2\}5\.28×10−25\.28\\times 10^\{\-2\}42\.00×10−12\.00\\times 10^\{\-1\}1\.88×10−11\.88\\times 10^\{\-1\}1\.16×10−21\.16\\times 10^\{\-2\}5\.81×10−25\.81\\times 10^\{\-2\}51\.67×10−11\.67\\times 10^\{\-1\}1\.53×10−11\.53\\times 10^\{\-1\}1\.34×10−21\.34\\times 10^\{\-2\}8\.02×10−28\.02\\times 10^\{\-2\}61\.43×10−11\.43\\times 10^\{\-1\}1\.32×10−11\.32\\times 10^\{\-1\}1\.07×10−21\.07\\times 10^\{\-2\}7\.52×10−27\.52\\times 10^\{\-2\}71\.25×10−11\.25\\times 10^\{\-1\}1\.17×10−11\.17\\times 10^\{\-1\}8\.06×10−38\.06\\times 10^\{\-3\}6\.45×10−26\.45\\times 10^\{\-2\}81\.11×10−11\.11\\times 10^\{\-1\}1\.04×10−11\.04\\times 10^\{\-1\}6\.97×10−36\.97\\times 10^\{\-3\}6\.27×10−26\.27\\times 10^\{\-2\}91\.00×10−11\.00\\times 10^\{\-1\}9\.29×10−29\.29\\times 10^\{\-2\}7\.11×10−37\.11\\times 10^\{\-3\}7\.11×10−27\.11\\times 10^\{\-2\}109\.09×10−29\.09\\times 10^\{\-2\}8\.54×10−28\.54\\times 10^\{\-2\}5\.48×10−35\.48\\times 10^\{\-3\}6\.03×10−26\.03\\times 10^\{\-2\}118\.33×10−28\.33\\times 10^\{\-2\}7\.79×10−27\.79\\times 10^\{\-2\}5\.42×10−35\.42\\times 10^\{\-3\}6\.51×10−26\.51\\times 10^\{\-2\}127\.69×10−27\.69\\times 10^\{\-2\}7\.22×10−27\.22\\times 10^\{\-2\}4\.69×10−34\.69\\times 10^\{\-3\}6\.09×10−26\.09\\times 10^\{\-2\}137\.14×10−27\.14\\times 10^\{\-2\}6\.72×10−26\.72\\times 10^\{\-2\}4\.27×10−34\.27\\times 10^\{\-3\}5\.98×10−25\.98\\times 10^\{\-2\}146\.67×10−26\.67\\times 10^\{\-2\}6\.23×10−26\.23\\times 10^\{\-2\}4\.34×10−34\.34\\times 10^\{\-3\}6\.51×10−26\.51\\times 10^\{\-2\}156\.25×10−26\.25\\times 10^\{\-2\}5\.87×10−25\.87\\times 10^\{\-2\}3\.77×10−33\.77\\times 10^\{\-3\}6\.03×10−26\.03\\times 10^\{\-2\}165\.88×10−25\.88\\times 10^\{\-2\}5\.58×10−25\.58\\times 10^\{\-2\}3\.06×10−33\.06\\times 10^\{\-3\}5\.20×10−25\.20\\times 10^\{\-2\}175\.56×10−25\.56\\times 10^\{\-2\}5\.32×10−25\.32\\times 10^\{\-2\}2\.31×10−32\.31\\times 10^\{\-3\}4\.16×10−24\.16\\times 10^\{\-2\}185\.26×10−25\.26\\times 10^\{\-2\}5\.07×10−25\.07\\times 10^\{\-2\}1\.91×10−31\.91\\times 10^\{\-3\}3\.64×10−23\.64\\times 10^\{\-2\}195\.00×10−25\.00\\times 10^\{\-2\}4\.82×10−24\.82\\times 10^\{\-2\}1\.81×10−31\.81\\times 10^\{\-3\}3\.62×10−23\.62\\times 10^\{\-2\}204\.76×10−24\.76\\times 10^\{\-2\}4\.60×10−24\.60\\times 10^\{\-2\}1\.66×10−31\.66\\times 10^\{\-3\}3\.49×10−23\.49\\times 10^\{\-2\}214\.55×10−24\.55\\times 10^\{\-2\}4\.42×10−24\.42\\times 10^\{\-2\}1\.22×10−31\.22\\times 10^\{\-3\}2\.67×10−22\.67\\times 10^\{\-2\}224\.35×10−24\.35\\times 10^\{\-2\}4\.26×10−24\.26\\times 10^\{\-2\}9\.28×10−49\.28\\times 10^\{\-4\}2\.13×10−22\.13\\times 10^\{\-2\}234\.17×10−24\.17\\times 10^\{\-2\}4\.06×10−24\.06\\times 10^\{\-2\}1\.11×10−31\.11\\times 10^\{\-3\}2\.67×10−22\.67\\times 10^\{\-2\}244\.00×10−24\.00\\times 10^\{\-2\}3\.89×10−23\.89\\times 10^\{\-2\}1\.07×10−31\.07\\times 10^\{\-3\}2\.68×10−22\.68\\times 10^\{\-2\}253\.85×10−23\.85\\times 10^\{\-2\}3\.74×10−23\.74\\times 10^\{\-2\}1\.07×10−31\.07\\times 10^\{\-3\}2\.78×10−22\.78\\times 10^\{\-2\}263\.70×10−23\.70\\times 10^\{\-2\}3\.60×10−23\.60\\times 10^\{\-2\}1\.00×10−31\.00\\times 10^\{\-3\}2\.71×10−22\.71\\times 10^\{\-2\}273\.57×10−23\.57\\times 10^\{\-2\}3\.51×10−23\.51\\times 10^\{\-2\}6\.14×10−46\.14\\times 10^\{\-4\}1\.72×10−21\.72\\times 10^\{\-2\}283\.45×10−23\.45\\times 10^\{\-2\}3\.34×10−23\.34\\times 10^\{\-2\}1\.11×10−31\.11\\times 10^\{\-3\}3\.23×10−23\.23\\times 10^\{\-2\}293\.33×10−23\.33\\times 10^\{\-2\}3\.23×10−23\.23\\times 10^\{\-2\}1\.03×10−31\.03\\times 10^\{\-3\}3\.08×10−23\.08\\times 10^\{\-2\}Second Moments02\.00×1002\.00\\times 10^\{0\}1\.91×1001\.91\\times 10^\{0\}8\.74×10−28\.74\\times 10^\{\-2\}4\.37×10−24\.37\\times 10^\{\-2\}15\.00×10−15\.00\\times 10^\{\-1\}4\.50×10−14\.50\\times 10^\{\-1\}4\.97×10−24\.97\\times 10^\{\-2\}9\.93×10−29\.93\\times 10^\{\-2\}22\.22×10−12\.22\\times 10^\{\-1\}2\.07×10−12\.07\\times 10^\{\-1\}1\.50×10−21\.50\\times 10^\{\-2\}6\.77×10−26\.77\\times 10^\{\-2\}31\.25×10−11\.25\\times 10^\{\-1\}1\.18×10−11\.18\\times 10^\{\-1\}6\.76×10−36\.76\\times 10^\{\-3\}5\.41×10−25\.41\\times 10^\{\-2\}48\.00×10−28\.00\\times 10^\{\-2\}7\.47×10−27\.47\\times 10^\{\-2\}5\.33×10−35\.33\\times 10^\{\-3\}6\.67×10−26\.67\\times 10^\{\-2\}55\.56×10−25\.56\\times 10^\{\-2\}4\.86×10−24\.86\\times 10^\{\-2\}6\.97×10−36\.97\\times 10^\{\-3\}1\.25×10−11\.25\\times 10^\{\-1\}64\.08×10−24\.08\\times 10^\{\-2\}3\.44×10−23\.44\\times 10^\{\-2\}6\.39×10−36\.39\\times 10^\{\-3\}1\.57×10−11\.57\\times 10^\{\-1\}73\.13×10−23\.13\\times 10^\{\-2\}1\.85×10−21\.85\\times 10^\{\-2\}1\.28×10−21\.28\\times 10^\{\-2\}4\.09×10−14\.09\\times 10^\{\-1\}82\.47×10−22\.47\\times 10^\{\-2\}2\.07×10−22\.07\\times 10^\{\-2\}4\.04×10−34\.04\\times 10^\{\-3\}1\.64×10−11\.64\\times 10^\{\-1\}92\.00×10−22\.00\\times 10^\{\-2\}1\.26×10−21\.26\\times 10^\{\-2\}7\.42×10−37\.42\\times 10^\{\-3\}3\.71×10−13\.71\\times 10^\{\-1\}101\.65×10−21\.65\\times 10^\{\-2\}9\.04×10−39\.04\\times 10^\{\-3\}7\.49×10−37\.49\\times 10^\{\-3\}4\.53×10−14\.53\\times 10^\{\-1\}111\.39×10−21\.39\\times 10^\{\-2\}8\.66×10−38\.66\\times 10^\{\-3\}5\.23×10−35\.23\\times 10^\{\-3\}3\.77×10−13\.77\\times 10^\{\-1\}121\.18×10−21\.18\\times 10^\{\-2\}−2\.12×10−3\-2\.12\\times 10^\{\-3\}1\.40×10−21\.40\\times 10^\{\-2\}1\.18×1001\.18\\times 10^\{0\}131\.02×10−21\.02\\times 10^\{\-2\}2\.56×10−32\.56\\times 10^\{\-3\}7\.65×10−37\.65\\times 10^\{\-3\}7\.49×10−17\.49\\times 10^\{\-1\}148\.89×10−38\.89\\times 10^\{\-3\}2\.26×10−32\.26\\times 10^\{\-3\}6\.63×10−36\.63\\times 10^\{\-3\}7\.45×10−17\.45\\times 10^\{\-1\}157\.81×10−37\.81\\times 10^\{\-3\}6\.45×10−36\.45\\times 10^\{\-3\}1\.37×10−31\.37\\times 10^\{\-3\}1\.75×10−11\.75\\times 10^\{\-1\}166\.92×10−36\.92\\times 10^\{\-3\}5\.53×10−35\.53\\times 10^\{\-3\}1\.39×10−31\.39\\times 10^\{\-3\}2\.01×10−12\.01\\times 10^\{\-1\}176\.17×10−36\.17\\times 10^\{\-3\}3\.53×10−33\.53\\times 10^\{\-3\}2\.64×10−32\.64\\times 10^\{\-3\}4\.28×10−14\.28\\times 10^\{\-1\}185\.54×10−35\.54\\times 10^\{\-3\}2\.21×10−32\.21\\times 10^\{\-3\}3\.33×10−33\.33\\times 10^\{\-3\}6\.01×10−16\.01\\times 10^\{\-1\}195\.00×10−35\.00\\times 10^\{\-3\}4\.34×10−34\.34\\times 10^\{\-3\}6\.56×10−46\.56\\times 10^\{\-4\}1\.31×10−11\.31\\times 10^\{\-1\}204\.54×10−34\.54\\times 10^\{\-3\}−3\.56×10−3\-3\.56\\times 10^\{\-3\}8\.09×10−38\.09\\times 10^\{\-3\}1\.78×1001\.78\\times 10^\{0\}214\.13×10−34\.13\\times 10^\{\-3\}−3\.49×10−4\-3\.49\\times 10^\{\-4\}4\.48×10−34\.48\\times 10^\{\-3\}1\.08×1001\.08\\times 10^\{0\}223\.78×10−33\.78\\times 10^\{\-3\}−6\.16×10−4\-6\.16\\times 10^\{\-4\}4\.40×10−34\.40\\times 10^\{\-3\}1\.16×1001\.16\\times 10^\{0\}233\.47×10−33\.47\\times 10^\{\-3\}9\.72×10−59\.72\\times 10^\{\-5\}3\.38×10−33\.38\\times 10^\{\-3\}9\.72×10−19\.72\\times 10^\{\-1\}243\.20×10−33\.20\\times 10^\{\-3\}−2\.92×10−3\-2\.92\\times 10^\{\-3\}6\.12×10−36\.12\\times 10^\{\-3\}1\.91×1001\.91\\times 10^\{0\}252\.96×10−32\.96\\times 10^\{\-3\}−2\.83×10−3\-2\.83\\times 10^\{\-3\}5\.79×10−35\.79\\times 10^\{\-3\}1\.96×1001\.96\\times 10^\{0\}262\.74×10−32\.74\\times 10^\{\-3\}−3\.85×10−3\-3\.85\\times 10^\{\-3\}6\.60×10−36\.60\\times 10^\{\-3\}2\.40×1002\.40\\times 10^\{0\}272\.55×10−32\.55\\times 10^\{\-3\}8\.05×10−58\.05\\times 10^\{\-5\}2\.47×10−32\.47\\times 10^\{\-3\}9\.68×10−19\.68\\times 10^\{\-1\}282\.38×10−32\.38\\times 10^\{\-3\}1\.71×10−31\.71\\times 10^\{\-3\}6\.69×10−46\.69\\times 10^\{\-4\}2\.81×10−12\.81\\times 10^\{\-1\}292\.22×10−32\.22\\times 10^\{\-3\}1\.58×10−31\.58\\times 10^\{\-3\}6\.43×10−46\.43\\times 10^\{\-4\}2\.89×10−12\.89\\times 10^\{\-1\}

## References

- \[1\]J\. Abate and W\. Whitt\(1992\)The fourier\-series method for inverting transforms of probability distributions\.Queueing systems10\(1\),pp\. 5–87\.Cited by:[§2](https://arxiv.org/html/2607.08091#S2.p1.5)\.
- \[2\]J\. Abate and W\. Whitt\(2006\)A unified framework for numerically inverting laplace transforms\.INFORMS Journal on Computing18\(4\),pp\. 408–421\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.p3.2)\.
- \[3\]B\. Ata, J\. M\. Harrison, and N\. Si\(2024\)Singular control of \(reflected\) brownian motion: a computational method suitable for queueing applications\.Queueing Systems108\(3\),pp\. 215–251\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[4\]B\. Ata, J\. M\. Harrison, and N\. Si\(2025\)Drift control of high\-dimensional reflected brownian motion: a computational method based on neural networks\.Stochastic Systems15\(2\),pp\. 111–146\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[5\]B\. Ata and Y\. Xu\(2025\)Dynamic control of stochastic matching systems in heavy traffic: an effective computational method for high\-dimensional problems\.arXiv preprint arXiv:2509\.00809\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[6\]B\. Ata and E\. Kaşıkaralar\(2025\)Dynamic scheduling of a multiclass queue in the halfin–whitt regime: a computational approach for high\-dimensional problems\.Management Science\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[7\]B\. Ata, W\. van Eekelen, and Y\. Zhong\(2025\)A computational method for solving the stochastic joint replenishment problem in high dimensions\.arXiv preprint arXiv:2511\.11830\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[8\]J\. Blanchet, X\. Chen, N\. Si, and P\. W\. Glynn\(2021\)Efficient steady\-state simulation of high\-dimensional stochastic networks\.Stochastic Systems11\(2\),pp\. 174–192\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.p3.2)\.
- \[9\]G\. Carleo and M\. Troyer\(2017\)Solving the quantum many\-body problem with artificial neural networks\.Science355\(6325\),pp\. 602–606\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p3.1)\.
- \[10\]G\. L\. Choudhury and D\. M\. Lucantoni\(1996\)Numerical computation of the moments of a probability distribution from its transform\.Operations Research44\(2\),pp\. 368–381\.Cited by:[§B\.1](https://arxiv.org/html/2607.08091#A2.SS1.p1.1)\.
- \[11\]J\. Dai and J\. M\. Harrison\(1992\)Reflected brownian motion in an orthant: numerical methods for steady\-state analysis\.The Annals of Applied Probability2\(1\),pp\. 65–86\.Cited by:[Appendix A](https://arxiv.org/html/2607.08091#A1.SS0.SSS0.Px4.p1.1),[§3](https://arxiv.org/html/2607.08091#S3.p2.4)\.
- \[12\]J\. Dai, M\. Miyazawa, and J\. Wu\(2014\)A multi\-dimensional srbm: geometric views of its product form stationary distribution\.Queueing Systems78\(4\),pp\. 313–335\.Cited by:[Appendix A](https://arxiv.org/html/2607.08091#A1.SS0.SSS0.Px4.p2.1),[Appendix A](https://arxiv.org/html/2607.08091#A1.SS0.SSS0.Px4.p2.5),[Appendix A](https://arxiv.org/html/2607.08091#A1.SS0.SSS0.Px4.p2.6),[§B\.2](https://arxiv.org/html/2607.08091#A2.SS2.p1.9),[§1](https://arxiv.org/html/2607.08091#S1.p2.14),[§1](https://arxiv.org/html/2607.08091#S1.p2.19),[§3](https://arxiv.org/html/2607.08091#S3.p2.4)\.
- \[13\]J\. Han, A\. Jentzen, and W\. E\(2018\)Solving high\-dimensional partial differential equations using deep learning\.Proceedings of the National Academy of Sciences115\(34\),pp\. 8505–8510\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[14\]J\. Han, A\. Jentzen,et al\.\(2017\)Deep learning\-based numerical methods for high\-dimensional parabolic partial differential equations and backward stochastic differential equations\.Communications in mathematics and statistics5\(4\),pp\. 349–380\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p2.1)\.
- \[15\]J\. M\. Harrison\(1978\)The diffusion approximation for tandem queues in heavy traffic\.Advances in Applied Probability10\(4\),pp\. 886–905\.Cited by:[§3](https://arxiv.org/html/2607.08091#S3.p2.4)\.
- \[16\]T\. mpmath development team\(2026\)Mpmath: a Python library for arbitrary\-precision floating\-point arithmetic \(version 1\.4\.0\)\.Note:http://mpmath\.org/Cited by:[§3](https://arxiv.org/html/2607.08091#S3.p2.4)\.
- \[17\]Y\. Qu, J\. Blanchet, and P\. Glynn\(2024\)Deep learning for computing convergence rates of markov chains\.Advances in Neural Information Processing Systems37,pp\. 84777–84798\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p1.1)\.
- \[18\]Y\. Qu, J\. Blanchet, and P\. Glynn\(2026\)Deep learning for markov chains: lyapunov functions, poisson’s equation, and stationary distributions\.Queueing Systems110\(1\),pp\. 10\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p1.1)\.
- \[19\]J\. Sirignano and K\. Spiliopoulos\(2018\)DGM: a deep learning algorithm for solving partial differential equations\.Journal of computational physics375,pp\. 1339–1364\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p3.1)\.
- \[20\]A\. Talbot\(1979\)The accurate numerical inversion of laplace transforms\.IMA Journal of Applied Mathematics23\(1\),pp\. 97–120\.Cited by:[§2](https://arxiv.org/html/2607.08091#S2.p1.5),[§3](https://arxiv.org/html/2607.08091#S3.p1.1)\.
- \[21\]J\. A\. C\. Weideman\(2006\)Optimizing talbot’s contours for the inversion of the laplace transform\.SIAM Journal on Numerical Analysis44\(6\),pp\. 2342–2362\.Cited by:[§2](https://arxiv.org/html/2607.08091#S2.p1.5)\.
- \[22\]E\. Weinan, J\. Han, and A\. Jentzen\(2021\)Algorithms for solving high dimensional pdes: from nonlinear monte carlo to machine learning\.Nonlinearity35\(1\),pp\. 278\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p3.1)\.
- \[23\]B\. Yuet al\.\(2018\)The deep ritz method: a deep learning\-based numerical algorithm for solving variational problems\.Communications in Mathematics and Statistics6\(1\),pp\. 1–12\.Cited by:[§1](https://arxiv.org/html/2607.08091#S1.SS0.SSS0.Px1.p3.1)\.