[Paper Review] MatDRAM: A pure-MATLAB Delayed-Rejection Adaptive Metropolis-Hastings Markov Chain Monte Carlo Sampler

Shashank Kumbhare, Amir Shahmoradi|arXiv (Cornell University)|Oct 8, 2020

Markov Chains and Monte Carlo Methods31 references4 citations

TL;DR

MatDRAM is a pure-MATLAB implementation of the Delayed-Rejection Adaptive Metropolis (DRAM) Markov Chain Monte Carlo (MCMC) algorithm designed for Bayesian inference, sampling, and integration in high-dimensional problems. It enables nearly fully automated, reproducible, and restartable MCMC simulations with adaptive proposal distributions and automatic monitoring of diminishing adaptation to ensure convergence.

ABSTRACT

Markov Chain Monte Carlo (MCMC) algorithms are widely used for stochastic optimization, sampling, and integration of mathematical objective functions, in particular, in the context of Bayesian inverse problems and parameter estimation. For decades, the algorithm of choice in MCMC simulations has been the Metropolis-Hastings (MH) algorithm. An advancement over the traditional MH-MCMC sampler is the Delayed-Rejection Adaptive Metropolis (DRAM). In this paper, we present MatDRAM, a stochastic optimization, sampling, and Monte Carlo integration toolbox in MATLAB which implements a variant of the DRAM algorithm for exploring the mathematical objective functions of arbitrary-dimensions, in particular, the posterior distributions of Bayesian models in data science, Machine Learning, and scientific inference. The design goals of MatDRAM include nearly-full automation of MCMC simulations, user-friendliness, fully-deterministic reproducibility, and the restart functionality of simulations. We also discuss the implementation details of a technique to automatically monitor and ensure the diminishing adaptation of the proposal distribution of the DRAM algorithm and a method of efficiently storing the resulting simulated Markov chains. The MatDRAM library is open-source, MIT-licensed, and permanently located and maintained as part of the ParaMonte library at https://github.com/cdslaborg/paramonte.

Motivation & Objective

• To develop a fully automated, user-friendly, and reproducible MCMC sampling toolbox in MATLAB for Bayesian inference and uncertainty quantification.
• To implement a variant of the Delayed-Rejection Adaptive Metropolis (DRAM) algorithm that ensures diminishing adaptation of the proposal distribution for convergence.
• To enable deterministic restart functionality for long-running MCMC simulations, especially on high-performance computing systems with time limits.
• To provide efficient storage and postprocessing of Markov chains with built-in convergence diagnostics and visualization tools.
• To ensure compatibility with the broader ParaMonte library for cross-language interoperability and shared postprocessing workflows.

Proposed method

• The sampler uses a multivariate normal (MVN) proposal distribution whose covariance is adaptively updated during simulation to improve sampling efficiency.
• It implements delayed rejection to allow multiple attempts at proposing new states, increasing the acceptance rate in complex posterior distributions.
• The algorithm monitors the Hellinger distance between successive proposal distributions to ensure diminishing adaptation, a key condition for ergodicity and convergence.
• An upper bound on the Total Variation Distance (TVD) between proposal distributions is computed and reported to assess adaptation progress.
• The library supports deterministic restarts by saving the state of the Markov chain and random number generator, ensuring bit-identical results across runs.
• Output chains are stored efficiently with metadata, and postprocessing tools are integrated for visualization and convergence analysis.

Experimental results

Research questions

• RQ1 How can a pure-MATLAB implementation of the DRAM algorithm be designed to ensure full reproducibility and restart capability in long MCMC simulations?
• RQ2 What techniques can be used to automatically monitor and enforce diminishing adaptation in adaptive MCMC samplers?
• RQ3 How effective is the Hellinger distance-based adaptation measure in detecting convergence and preventing non-ergodic behavior?
• RQ4 Can a fully automated, user-friendly MCMC toolbox be built in MATLAB that supports high-dimensional Bayesian inference with minimal user intervention?
• RQ5 How does the integration of MatDRAM with the ParaMonte library enable cross-language compatibility and unified postprocessing of MCMC results?

Key findings

• The MatDRAM library successfully implements a fully automated, deterministic, and restartable DRAM sampler in pure MATLAB, ensuring bit-identical results across interrupted and resumed simulations.
• The adaptation measure based on Hellinger distance provides a reliable, computationally efficient upper bound on Total Variation Distance (TVD), which decays progressively during simulation, indicating successful diminishing adaptation.
• The method of monitoring adaptation via Hellinger distance enables automatic detection of non-convergent chains, where adaptation fails to diminish.
• The library supports high-dimensional sampling (evidenced by 4D MVN examples) and efficiently stores Markov chains with full metadata for postprocessing.
• The integration with the ParaMonte library allows MatDRAM to interoperate with C++, Fortran, and Python-based samplers, enabling unified analysis of MCMC output.
• The tool is open-source, MIT-licensed, and hosted at https://github.com/cdslaborg/paramonte, with public access via MATLAB FileExchange and comprehensive documentation.