Introduction to LangevinFlow

Overview

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

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:

  1. Run a short pilot chain.
  2. If the acceptance rate is too high (>0.7), the steps are too small — increase step_size.
  3. If it is too low (<0.4), the steps are too large — decrease step_size.
  4. 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