Optimal Thinning of MCMC Output

Abstract

The use of heuristics to assess the convergence and compress the output of Markov chain Monte Carlo can be sub-optimal in terms of the empirical approximations that are produced. Typically a number of the initial states are attributed to 'burn in' and removed, while the remainder of the chain is 'thinned' if compression is also required. In this paper, we consider the problem of retrospectively selecting a subset of states, of fixed cardinality, from the sample path such that the approximation provided by their empirical distribution is close to optimal. A novel method is proposed, based on greedy minimisation of a kernel Stein discrepancy, that is suitable when the gradient of the log-target can be evaluated and approximation using a small number of states is required. Theoretical results guarantee consistency of the method and its effectiveness is demonstrated in the challenging context of parameter inference for ordinary differential equations. Software is available in the Stein Thinning package in Python, R and MATLAB.

Gradient-Free Kernel Stein Discrepancy

“…Proof. This is Lemma 4 in Riabiz et al [2022], specialised to the case where samples are independent and identically distributed. Although not identical to the statement in Riabiz et al [2022], one obtains this result by following an identical argument and noting that the expectation of k q (x i , x j ) is identically 0 when i = j (due to independence of x i and x j ), so that bounds on these terms are not required.…”

“…The proof builds on earlier work in Riabiz et al [2022]. Note that the optimal weights w * can be computed without the normalisation constant of p, by solving a constrained quadratic programme at cost O(n 3 ).…”

“…The direct application of Stein importance sampling to this task requires the numerical calculation of sensitivities of the differential equation at each of the n samples that are to be re-weighted. Aside from simple cases where automatic differentiation or adjoint methods can be used, the stable computation of sensitivities can form a major computational bottleneck; see [Riabiz et al, 2022]. In contrast, our approach required a fixed number of gradient computations to construct a Laplace approximation 3 , independent of the number n of samples required; see Appendix B.3 for detail.…”

“…This appears to be a challenging problem but we hope to address it in future work. In addition, although we focused on Stein importance sampling, our methodology offers the possibility to construct gradient-free versions of other related algorithms, including the Stein points algorithms of Chen et al [2018Chen et al [ , 2019, and the Stein thinning algorithm of Riabiz et al [2022].…”

“…To date, test functions have been constructed by combining derivatives of the statistical model with reproducing kernels [Chwialkowski et al, 2016, Gorham and Mackey, 2017, Gong et al, 2021a, random features [Huggins and Mackey, 2018], diffusion coefficients and functions with bounded derivatives , neural networks [Grathwohl et al, 2020], and polynomials [Chopin and Ducrocq, 2021]. The resulting discrepancies have been shown to be powerful statistical tools, with diverse applications including parameter inference , Matsubara et al, 2022, goodness-of-fit testing [Jitkrittum et al, 2017, Fernandez et al, 2020, and sampling [Liu and Lee, 2017, Chen et al, 2018, Riabiz et al, 2022, Hodgkinson et al, 2020, Fisher et al, 2021. However, one of the main drawbacks of these existing works is the requirement that derivatives both exist and can be computed.…”