The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction to LangevinFlow

Overview

LangevinFlow implements two Markov chain Monte Carlo (MCMC) samplers based on overdamped Langevin diffusion:

ula() — the Unadjusted Langevin Algorithm, an Euler-Maruyama discretization of the SDE. Cheap per step, but introduces a discretization bias that vanishes only as the step size goes to zero.
mala() — the Metropolis-Adjusted Langevin Algorithm, which corrects that bias with an accept/reject step. Asymptotically exact.

Both target \(\pi(x) \propto \exp(-\beta U(x))\) for a user-supplied potential \(U\). The user passes the gradient \(\nabla U\) (and, for MALA, \(U\) itself).

Sign convention

The Langevin SDE is

\[ dX_t = -\nabla U(X_t)\,dt + \sqrt{2/\beta}\,dW_t, \]

so the user supplies \(U = -\log \pi\) (the potential / negative log density). Get this wrong and the chain runs away from the mode rather than toward it — so it’s worth a sanity check on a simple target before trusting any real run.

Example 1: standard Gaussian

library(LangevinFlow)

# Target: N(0, I_2). U(x) = 0.5 * ||x||^2.
U      <- function(x) 0.5 * sum(x^2)
grad_u <- function(x) x

set.seed(1)
fit_ula <- ula(init_x = c(3, -3), grad_u = grad_u,
               step_size = 0.05, n_iter = 5000, burn_in = 1000)
summary(fit_ula)
#> <langevin_chain>
#>   Algorithm:       ULA
#>   Dimension:       2
#>   Iterations:      5000 (burn-in: 1000)
#>   Step size:       0.05
#>   Inverse temp.:   1
#>   Elapsed:         0.005 sec
#> 
#> Per-coordinate summaries:
#>        mean     sd   2.5%      50% 97.5%
#> x1 -0.13660 1.0040 -2.109 -0.14610 1.815
#> x2  0.07887 0.9981 -1.892  0.08295 2.077

plot(fit_ula)

ULA trace plots

The same target with MALA:

set.seed(1)
fit_mala <- mala(init_x = c(3, -3), U = U, grad_u = grad_u,
                 step_size = 0.4, n_iter = 5000, burn_in = 1000)
fit_mala$acceptance_rate
#> [1] 0.9068
summary(fit_mala)
#> <langevin_chain>
#>   Algorithm:       MALA
#>   Dimension:       2
#>   Iterations:      5000 (burn-in: 1000)
#>   Step size:       0.4
#>   Inverse temp.:   1
#>   Acceptance rate: 0.907
#>   Elapsed:         0.010 sec
#> 
#> Per-coordinate summaries:
#>         mean     sd   2.5%       50% 97.5%
#> x1 0.0002337 0.9776 -1.946 -0.008739 1.911
#> x2 0.0346300 0.9726 -1.844  0.027750 1.990

Example 2: a correlated Gaussian

Sigma     <- matrix(c(1, 0.8, 0.8, 1), 2, 2)
Sigma_inv <- solve(Sigma)
mu_true   <- c(2, -1)

U      <- function(x) 0.5 * as.numeric(t(x - mu_true) %*% Sigma_inv %*% (x - mu_true))
grad_u <- function(x) as.numeric(Sigma_inv %*% (x - mu_true))

set.seed(2)
fit <- mala(c(0, 0), U, grad_u, step_size = 0.4,
            n_iter = 8000, burn_in = 2000)
colMeans(fit$samples)
#> [1]  1.948857 -1.059934
cov(fit$samples)
#>           [,1]      [,2]
#> [1,] 1.0760520 0.8517098
#> [2,] 0.8517098 1.0478167

Choosing the step size

For MALA in \(d\) dimensions, the optimal acceptance rate is known to be roughly 0.574 as \(d \to \infty\) (Roberts & Rosenthal, 1998). A practical recipe:

Run a short pilot chain.
If the acceptance rate is too high (>0.7), the steps are too small — increase step_size.
If it is too low (<0.4), the steps are too large — decrease step_size.
Iterate until it lands in the 0.5–0.6 range.

For ULA there is no accept/reject feedback, so step-size choice depends on the curvature of \(U\): too large and the chain diverges (the package raises an error); too small and mixing is slow.

References

Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4), 341-363.
Roberts, G. O., & Rosenthal, J. S. (1998). Optimal scaling of discrete approximations to Langevin diffusions. JRSS-B, 60(1), 255-268.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.