In mathematics, the Leimkuhler-Matthews method (or LM method in its original paper ) is an algorithm for finding discretized solutions to the Brownian dynamics
d X = − ∇ V ( X ) d t + σ d W , {\displaystyle \mathrm {d} X=-\nabla V(X)\,\mathrm {d} t+\sigma \,\mathrm {d} W,}
where σ > 0 {\displaystyle \sigma >0} is a constant, V ( X ) {\displaystyle V(X)} is an energy function and W ( t ) {\displaystyle W(t)} is a Wiener process. This stochastic differential equation has solutions (denoted X ( t ) ∈ R N {\displaystyle X(t)\in \mathbb {R} ^{N}} at time t {\displaystyle t} ) distributed according to π ( X ) ∝ exp ( − V ( x ) ) {\displaystyle \pi (X)\propto \exp(-V(x))} in the limit of large-time, making solving these dynamics relevant in sampling-focused applications such as classical molecular dynamics and machine learning. Given a time step Δ t > 0 {\displaystyle \Delta t>0} , the Leimkuhler-Matthews update scheme is compactly written as
X t + Δ t = X t − ∇ V ( X t ) Δ t + σ Δ t 2 ( R t + R t + Δ t ) , {\displaystyle X_{t+\Delta t}=X_{t}-\nabla V(X_{t})\Delta t+\sigma {\frac {\sqrt {\Delta t}}{2}}\,(R_{t}+R_{t+\Delta t}),}
with initial condition X 0 := X ( 0 ) {\displaystyle X_{0}:=X(0)} , and where X t ≈ X ( t ) {\displaystyle X_{t}\approx X(t)} . The vector R t {\displaystyle R_{t}} is a vector of independent normal random numbers redrawn at each step so E [ R t ⋅ R s ] = N δ t s {\displaystyle {\text{E}}[R_{t}\cdot R_{s}]=N\delta _{ts}} (where E [ ∙ ] {\displaystyle {\text{E}}[\bullet ]} denotes expectation). Despite being of equal cost to the Euler-Maruyama scheme (in terms of the number of evaluations of the function ∇ V ( X ) {\displaystyle \nabla V(X)} per update), given some assumptions on Δ t , V ( X ) {\displaystyle \Delta t,\,V(X)} and f ( X ) {\displaystyle f(X)} solutions have been shown to have a superconvergence property
| E [ f ( X t ) ] − E [ f ( X ( t ) ) ] | ≤ C 1 e − λ t Δ t + C 2 Δ t 2 {\displaystyle |{\text{E}}[f(X_{t})]-{\text{E}}[f(X(t))]|\leq C_{1}e^{-\lambda t}\Delta t+C_{2}\Delta t^{2}}
for constants C k ≥ 0 , λ > 0 {\displaystyle C_{k}\geq 0,\,\lambda >0} not depending on t {\displaystyle t} . This means that as t {\displaystyle t} gets large we obtain an effective second order with Δ t 2 {\displaystyle \Delta t^{2}} error in computed expectations. For small time step Δ t {\displaystyle \Delta t} this can give significant improvements over the Euler-Maruyama scheme, at no extra cost.
Discussion
Comparison to other schemes The obvious method for comparison is the Euler-Maruyama scheme as it has the same cost, requiring one evaluation of ∇ V ( X ) {\displaystyle \nabla V(X)} per step. Its update is of the form
… excerpt ends here. Continue reading the full article.

