Home People Software Projects Alumni Links

Computational Stochastics

Our group studies Monte Carlo methods computational methods for random systems. The common theme in these projects is the use of mathematical analysis to understand and design new Monte Carlo methods.

Time stepping methods for stochastic differential equations

A stochastic differential equation may be written as

A time stepping method is a way to generate sample paths that can be used for Monte Carlo. The best known methods are the (forward) Euler method and Milstein's method. Milstein's method is more accurate than the Euler method in a strong pathwise sense but not in a weaker statistical sense more appropriate for Monte Carlo. We have been exploring several related issues -- finding alternative error measures that are statistical but still measure path properties, and finding alternative time stepping methods that have the accuracy of Milstein's method but not some of its drawbacks.

In particular, we found a two stage Runge Kutta method that achieves the accuracy of Milstein's method without using analytical derivatives of the coefficients and without needing to sample the Levi area integral (a random variable defined by a stochastic integral of two independent Brownian motions). This Runge Kutta method would not be possible without our new error measure. We are currently working on a different error measure that is closer to the traditional strong error. This is the expected difference between two paths, in an optimal coupling. The so called ``coupling distance'' is related to the Wasserstein, or Monge- Kantorovich-Rubinstein metric.

Computing slow relaxation modes

Many complex stochastic dynamical systems have subtle and very slow relaxation modes. Dynamical glasses and the glass transition are one important class of examples. We are developing a method that more accurately computes the second eigenvalue of the operator, which governs the rate at which the Markov chain statistics converges to equilibrium. So far our method is much more accurate and less noisy than Prony's method, the standard approach for these problems.

For technical details see: K. Gade, Numerical Estimation of the Second Largest Eigenvalue of Reversible Markov Transition Matrix , PhD Thesis, NYU, January, 2009.

Dynamic coupling control variates

Control variates are a way to reduce the variance in a Monte Carlo computation by subtracting out the known answer to a simpler problem. We were interested in model one dimensional lattice systems from statistical mechanics that would satisfy detailed balance in equilibrium. Out of equilibrium, there is net energy flux through the chain that we want to estimate by Monte Carlo simulation. The problem is that long chains with a fixed temperature difference at the ends locally are close to equilibrium. This makes is hard to estimate the relatively small energy flux in the presence of statistical noise.

Our method attempts to reduce the variance by coupling the dynamic simulation to another simulation that preserves an approximation to the local equilibrium distribution. The coupling is achieved using a Metropolis style rejection strategy. Computational experiments show improvements in accuracy by roughly a factor of two, which is less than should be possible with such couplings. We hope that a more sophisticated coupling to a more accurate approximate distribution will give much better accuracy.

For technical details see: J. Goodman and K. Lin, Dynamic control variates for near equilibrium systems, preprint.

Rare events without large deviations

Sometimes it is necessary to estimate the probability of very unlikely events. For example, we may want to know the failure probability of a very reliable system. Chemical reactions are rare on the time scale of atomic motion, but do happen on a macroscopic time scale. In the case of a reliable system or a slow chemical reaction, direct simulation is not an effective way to find the rare events because most of the simulation time does not see them.

Rare event simulation refers to Monte Carlo methods that find rare events more efficiently than direct simulation. Systematic methods for rare event simulation often are based on the mathematical theory of large deviations (for which Courant Institute mathematician S.R.S Varadhan received the Abel Prize). Large deviation theory essentially asks the user to find the most likely mechanism for producing rare events analytically before starting the computation.

Our group is seeking ways to find rare events without solving the large deviation problem first. The picture below was produced using a histogram bisection strategy. Suppose X represents a configuration of some complex random system and we want to know the probability that f(X) > A, when this event is very small. Our method proceeds by repeated bifurcation of the histogram. The first step finds the median of the histogram, which is the number M1 so that the probability of f(X)>M1 is equal to .5. The next stage is to refine this top half to find M2 so that the probability that f(X)>m2 is .25. We do this not by direct simulation, but by a Metropolis type sampler that rejects configurations X that fall outside f(X)>.5. In the figure below, this bifurcation process was repeated 20 times, giving paths whose probability is less than one in a million. The whole process took well under a million path simulations. The large deviation theory for this problem is known and the path below is typical of what large deviations predicts such paths should look like.