BAOAB underdamped Langevin splitting

Use the equation and conventions of Underdamped Langevin equation and canonical Gibbs density. Given a valid state $(q_n,p_n)\in {\mathbb{R}}^{2d}$, a step $h>0$, and a fresh independent vector ${\xi }_n\sim N_d(0,I_d)$, put $a=e^{-\gamma h}$ and execute exactly these five substeps in order:

1. B half kick: $p^{(1)}=p_n-\frac{h}{2}\mathrm{\nabla }V(q_n)$.
2. A half drift: $q^{(1)}=q_n+\frac{h}{2}{M}^{-1}{p}^{(1)}$.
3. O full thermostat: $p^{(2)}=a{p}^{(1)}+\sqrt{{\beta }^{-1}(1-a^2)}{M}^{1\mathrm{/}2}{\xi }_n$.
4. A half drift: $q_{n+1}=q^{(1)}+\frac{h}{2}{M}^{-1}{p}^{(2)}$.
5. B half kick: $p_{n+1}=p^{(2)}-\frac{h}{2}\mathrm{\nabla }V(q_{n+1})$.

The output is $(q_{n+1},p_{n+1})$. If any required gradient is undefined or any input violates $h>0$, $\beta >0$, $\gamma >0$, or positive definiteness of $M$, return failure without an updated state. For a differentiable $V:{\mathbb{R}}^d\to \mathbb{R}$ and finite valid state, both gradient calls are defined; reusing the stored value $\mathrm{\nabla }V(q_{n+1})$ as the next step's first force means one new force evaluation per subsequent step. The random-vector cost is one fresh $d$-dimensional standard normal per step.

The O step is the exact momentum Ornstein–Uhlenbeck transition proved in Exact Ornstein–Uhlenbeck momentum transition and covariance. When each ${\xi }_n$ is independent of the current state and all earlier draws, the conditional law of the output given the current state depends only on that state; this is the Markov transition defined by the algorithm. Its procedure definition makes no claim of general invariant-law existence, uniqueness, or bias order.

A probability ${\pi }_h$ is invariant for this transition when one update from a state with law ${\pi }_h$, using an independent fresh normal vector, again has law ${\pi }_h$. For such an invariant probability ${\pi }_h$ and target Gibbs probability $\pi$, invariant-distribution bias means the difference between expectations under ${\pi }_h$ and $\pi$ for a specified observable. Burn-in error is the difference between the law at a finite simulation time and ${\pi }_h$. Correlation error is the effect of dependence among successive outputs on the variance of a sample average, and sampling error is the random deviation of that sample average from its expectation under the simulation law. These are separate quantities: the present method gives no general bound for any of them. The Gaussian draws are mathematical exact-normal inputs. Finite-bit random-number generation error means the discrepancy between the implemented finite-precision draw law and the exact normal law specified here; it is a separate implementation error and is not covered by this algorithm.