NUTS for Truncated Hamiltonian Monte Carlo

Background

This project is about sampling from high dimensional log-concave distributions with Hamiltonian Monte Carlo (HMC) algorithm. In particular, the student will adjust and implement the dynamic stopping criterion of No-U-Turn Sampler (NUTS) [1] to the case of truncated by a convex polytope sampling. In the sequel, she/he will use numerical approximation of the gradient of -log[p(x)] in HMC step, where p(x) is the PDF of the density we wish to sample from. The aim is to provide an efficient zero-order HMC sampler, that will not require an evaluation oracle for the gradient by the user. The student's code will be based on volesti's implementation of HMC algorithm. NUTS criterion is used to improve the mixing rate of HMC and to automate the parameterization of the random walk.

Related work

Only volesti package provides an open source implementation for sampling from multivariate distributions truncated by a polytope with HMC algorithm. It employs several integrators for the ODE that HMC step requires. The current implementation does not use the NUTS criterion neither for truncated nor for the non-truncated case. The NUTS criterion is used also for non-truncate distributions by stan, a state-of-the-art platform for statistical modelling and high-performance statistical computation.

An alternative to leapfrog set of integrators is presented in this paper and could be extended in the NUTS case.

Details of your coding project

The student (a) will implement the NUTS criterion for HMC algorithm in C++, (b) will implement several gradient approximation methods, (c) evaluate and tune (a), (b) for several ODE solvers that volesti provides and (d) update volesti's documentation for the new implementations. The proposed programming language is C++ because the implemented code will have to be added to volesti package and GeomScale project in general. Of course, the implementation will be based on the current software of volesti, so the basic geometrical concepts (polytopes, boundary reflections, ODE solvers) are ready to be used.

Expected impact

This project is expected to be very important to GeomScale project as it is about efficient sampling from log-concave or general distributions which appear in many applications. Moreover, it will be an important stepping stone towards high dimensional MCMC sampling with open source software, addressing some important applications in biology and bayesian inference.

[1] The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Matthew D. Hoffman, Andrew Gelman

Mentors

Students, please contact both mentors below after completing at least one of the tests below.

• Vissarion Fisikopoulos <vissarion.fisikopoulos at gmail.com> is an expert in mathematical software, computational geometry and optimization, and has previous GSOC mentoring experience with Boost C++ libraries (2016-2020) and the R-project (2017-2020).

• Marios Papachristou < papachristoumarios at gmail.com > is a PhD student in the Computer Science Department at Cornell University. His primary research interests lie within the field of Data Science. He has previous experience in GSoC 2018 and 2020 as a student under Org. FOSS and GeomScale. He was GSoC mentor in GSoC 2019.

• Elias Tsigaridas <elias.tsigaridas at inria.fr> is an expert in computational nonlinear algebra and geometry with experience in mathematical software. He has contributed to the implementation, in C and C++, of several solving algorithms for various open source computer algebra libraries and has previous GSOC mentoring experience with the R-project (2019) and Geomscale (2020).

Tests

Students, please do one or more of the following tests before contacting the mentors above.

• Easy: Download, compile and run a simple sampling example with both C++ and R interfaces of volesti. For example, you can sample points from a 100-dimensional cube using the HMC algorithm implemented in volesti for various distributions.

• Medium: Add counters in C++ code of volesti to count the average number of reflections per HMC step. For different (a) integration time, (b) step length in ODE solver, (c) dimension, compute the PSRF/run-time and Effective Sample Size (ESS)/run-time for a large number of samples and report.

• Hard: Implement in C++ an approximate and an exact oracle for the gradient of the Spherical Gaussian distribution with unit variance and repeat the experiments requested by the Medium test. Compare the performances of two cases and report.