Algorithms for integration of stochastic differential equations using parallel optimized sampling in the Stratonovich calculus

Abstract

A variance reduction method for stochastic integration of Fokker–Planck equations is derived. This unifies the cumulant hierarchy and stochastic equation approaches to obtaining moments, giving a performance superior to either. We show that the brute force method of reducing sampling error by just using more trajectories in a sampled stochastic equation is not the best approach. The alternative of using a hierarchy of moment equations is also not optimal, as it may converge to erroneous answers. Instead, through Bayesian conditioning of the stochastic noise on the requirement that moment equations are satisfied, we obtain improved results with reduced sampling errors for a given number of stochastic trajectories. The method used here converges faster in time-step than Ito–Euler algorithms. This parallel optimized sampling (POS) algorithm is illustrated by several examples, including a bistable nonlinear oscillator case where moment hierarchies fail to converge.

Introduction

Fokker–Planck equations or FPEs  [1], [2], combining first order drift and second order diffusion terms, are among the most universally encountered equations in physics. First obtained in work of Fokker, Planck and Langevin [3] on particle motion in a random environment, their use has now extended to many areas of physics and other sciences [4], [5], [6]. They are equivalent to stochastic differential equations or SDEs [7], [8], [9]. Functional Fokker–Planck equations are equivalent to the increasingly useful stochastic partial differential equations [10].

In this work, we derive and analyze a novel, hybrid approach to solving these equations. Already, there are a variety of methods [5]. The most common, due to its ease of use, is the direct solution of the equivalent stochastic equation through numerical simulation. Yet stochastic methods are greatly restricted by sampling error. As a result, the random error is no better than $1/\sqrt{N}$ for $N$ trajectories. While there are methods to reduce the error due to time discretization [11], [12], [13], [14], [15], variance reduction to reduce sampling errors remains problematic [16], [17].

Another, apparently more efficient method, is to use a hierarchy of moment equations [18], [19], [20]. This approach requires truncation, as the hierarchy is generically infinite. The truncation approach has a severe drawback, however: truncated hierarchies of moment equations are fundamentally inconsistent, and can have multiple solutions [21]. As a result, moment equation hierarchies often may converge to an incorrect answer  [22]. This is a very severe limitation. A truncated hierarchy could appear to converge, yet the answer obtained may be wrong.

We introduce an alternative approach for solving Fokker–Planck equations called parallel optimized sampling (POS). This unifies the two earlier approaches, giving the efficiency of moment equations while retaining the correctness of stochastic methods. In the present paper, we show how to obtain a POS algorithm with moment equations from the Stratonovich form of a stochastic equation  [8]. This is a common form of equation for a physical stochastic process.

Stratonovich equations have the merit that they are well-suited to numerical integration. They are the broad-band limit of a colored noise–i.e., finite band-width–stochastic process. They follow the rules of standard calculus, and algorithms similar to the usual ordinary differential equation methods can be used. In principle, one can transform these to the Ito form, which is preferred by mathematicians, and then use methods designed for Ito equations. This approach is treated here  [23]. However, there can be numerical advantages to using algorithms that are close to ordinary methods of calculus.

In order to demonstrate how these Stratonovich-POS methods work, we analyze the application of this method to some well-known problems: the linear oscillator, the Kubo oscillator, stochastic motion in a bistable well, and a complex-valued laser equation. In the bistable case, the truncated moment hierarchy method gives incorrect moments even when taken to very high order. This illustrates the advantage of the hybrid POS approach. The complex-valued example demonstrates that these methods are applicable in more than one dimension.

We show that POS methods, unifying moment and stochastic methods, not only converge to the correct result, but give orders of magnitude greater accuracy than stochastic methods with the same trajectory number. Rather than trying to obtain precision with more independent trajectories, it is greatly advantageous to use the POS approach with fewer trajectories.

Sampled Monte Carlo methods  [24] are the mathematical foundation of many other methods in statistical physics and quantum many-body theory. Hierarchies of moments are also the basis of ‘dynamical mean-field’  [25], ‘cluster expansion’  [26] or ‘multi-configurational [27] types of approach. Accordingly, it is important to understand such methods. These approaches appear to have the same fundamental limitation, i.e. that direct simulations have large sampling errors, while nonlinear hierarchies of equations may have multiple, incorrect solutions. Thus, while we focus on stochastic equations here, it is not impossible that POS techniques could be useful for treating other statistical problems.