Fixed-Horizon Error Rates for Array-RQMC Metropolis-Hastings

Introduction

Monte Carlo (MC) is a standard approach for approximating target expectations, which we denote by $\Pi(f)=\int_{\mathcal X}f(\bm{x})\pi(\bm{x})\,\,\mathrm d\bm{x}$, where $\pi$ is a target density on a state space $\mathcal X$ and $f$ is a test function. In many applications, particularly in Bayesian computation, direct sampling from the target distribution $\Pi$ is infeasible and $\pi$ is known only up to a normalizing constant. Markov chain Monte Carlo (MCMC), and in particular the Metropolis–Hastings (MH) algorithm, addresses this problem by constructing a Markov chain with invariant distribution $\Pi$ and estimating $\Pi(f)$ from its sample path; see, for example, [1]. In standard Monte Carlo integration, the driving inputs are independent and identically distributed (IID), whereas quasi-Monte Carlo (QMC) replaces them by more evenly distributed deterministic point sets with the aim of obtaining more accurate numerical estimates. Randomized QMC (RQMC) randomizes such point sets while preserving their uniformity properties and permits error assessment through independent replications [2]. This naturally raises the question of whether QMC uniformity can also improve the accuracy of practical MH algorithms.

One approach is Markov chain quasi-Monte Carlo (MCQMC), initiated by [3], which replaces the successive IID inputs along one chain trajectory by a completely uniformly distributed (CUD) or weakly CUD sequence. Numerical studies [5, 4] have reported substantial variance reductions and, for several Gibbs samplers, empirical convergence rates faster than Monte Carlo. For MH, however, the reported gains are often more modest. The difficulty is that low-discrepancy structure must survive a long nonlinear composition with an acceptance discontinuity at every transition. Available MCQMC theory establishes consistency under suitable conditions [6, 7], while convergence-rate results rely on existence arguments [8, 9]; in particular, they do not provide rates for the explicit CUD constructions commonly used in practice.

Array-RQMC instead imposes QMC uniformity across a particle population at each fixed transition. The randomized construction was introduced in [10]; mapping alternatives and convergence behavior were studied further in [11], and the method, sorting strategies, and known convergence results were reviewed in [12]. It has since been applied in a range of simulation settings, including option pricing, stochastic reaction networks, and, more recently, walk-on-spheres simulation [13, 15, 14]. It evolves $N$ particles in parallel: at each time, it sorts the population, matches the particles to a freshly randomized low-discrepancy point set, and advances them simultaneously. Its fresh randomization at each transition prevents the nominal RQMC dimension from growing with $T$. For MH, this population structure, which is local in time, turns the acceptance discontinuity into an indicator integration problem of fixed dimension at each step. Existing convergence theory provides two useful benchmarks. [16] established an $o(N^{-1})$ variance result for randomized sequential QMC under regularity conditions, with corresponding Array-RQMC conclusions obtained through Hilbert ordering of multidimensional states. Earlier, [10] proved the sharper explicit rate $\mathcal O(N^{-3/2})$ for one-dimensional monotone inversion transitions under a square-stratified sampling construction and additional regularity assumptions. These results do not directly provide a rate for the usual proposal-and-accept/reject implementation of MH, which we analyze through the geometry of local indicator events.

Because Array-RQMC is intrinsically a population method, its natural MC benchmark is $N$ independent copies of the same chain evolved over the same finite horizon. For a fixed horizon $T$, averaging these copies has MSE $\mathcal{O}_T(N^{-1})$ for the corresponding expectation under a finite second moment. The question addressed here is whether the dependence deliberately introduced by Array-RQMC can improve this rate in $N$ for a standard MH transition.

Our analysis concerns sampling error relative to the finite-horizon expectation $\mu_T(f)$ with $T$ fixed. The Array-RQMC estimator is unbiased for this expectation; unless the chain starts in stationarity, the estimation of $\Pi(f)$ may contain an additional bias independent of $N$. In one dimension, explicit variation conditions give $\mathcal O_T(N^{-4/3})$ for bounded-variation test functions. For product targets, coordinatewise Array-RQMC Metropolis-within-Gibbs (MWG) retains this exponent for additive marginal bounded-variation tests. This is a dimension-independent exponent, but not a bound with dimension-uniform constants. With Gaussian proposals, an extension to estimating the expectation of polynomials has the additional factor $(\log N)^\ell$ for maximum polynomial degree $\ell$. For joint $p$-dimensional MH, suitable Hardy–Krause variation conditions and uniform boundary-tube bounds yield a conditional $\mathcal O_T(N^{-1-1/(p+2)})$ bound. These hypotheses do not follow from Hilbert sorting and are not verified for the joint examples.

Numerical experiments support the verified rates and contrast this product construction with joint Hilbert-sorted MH in increasing dimension.

§ 2 reviews Array-RQMC and Hilbert ordering; § 3 gives the MH algorithm, the multidimensional geometric result and the verifiable one-dimensional and product-target theory. § 4 reports experiments, and § 5 discusses the remaining questions. The appendices prove the multidimensional and one-dimensional fixed-horizon theorems. The supplementary materials follow the one-dimensional argument: supporting lemmas, verification tools and target classes, the product-target MWG reduction and polynomial-test extension, the scrambled-net inversion extension and horizon dependence. MATLAB scripts (MATLAB R2024b) are available from https://github.com/Jiarui-Du/Array-RQMC-MH.

Background and Array-RQMC

Notation

For a joint MH update in state dimension $p$, we consider proposals generated from $p$ uniform variables. $D$ denotes the dimension of each Array-RQMC point set: $D=p+2$ for joint MH, including the acceptance and sorting/matching coordinates, and $D=3$ for each coordinatewise MWG update in § 3.6. Vectors and vector-valued maps are written in bold, for example $\bm{x}=(x_1,\ldots,x_p)$; scalar components and one-dimensional state variables are written in ordinary type. We write $\lambda_D$ for the $D$-dimensional Lebesgue measure and $\mathbf 1_E$ for the indicator of a set $E$. For a real random variable $Z$, write $\|Z\|_2:=(\mathbb E|Z|^2)^{1/2}$, and for a sigma-field $\mathcal G$, write $\|Z\|_{2\mid\mathcal G}:=(\mathbb E[|Z|^2\mid\mathcal G])^{1/2}$. In our results this expectation is $\mathbb E_{\mathrm{RQMC}}$, taken over the point-set randomizations (§ 3.2). For a deterministic function $g$ on $[0,1)^D$, $\|g\|_2$ denotes its $L^2(\lambda_D)$ norm.

For a real-valued function $h$ on an interval $\mathcal X \subseteq\mathbb R$, let $V_{\mathcal X}(h):= \sup_{\substack{r\ge1\\x_0<\cdots<x_r,\;x_\ell\in\mathcal X}} \sum_{\ell=1}^r |h(x_\ell)-h(x_{\ell-1})|,$ and write ${\rm BV}(\mathcal X)$ for the class of bounded variation functions satisfying $V_{\mathcal X}(h)<\infty$. Constants hidden in $\mathcal{O}_T(\cdot)$ depend on the fixed horizon $T$, but not on $N$.

Array-RQMC and Hilbert ordering

Array-RQMC evolves $N$ copies of a Markov chain $\bm{X}_j=\bm{\varphi}_j(\bm{X}_{j-1},\bm{\Theta}_j)$ using a fresh randomized low-discrepancy point set at each transition [10, 11, 12], where $\bm{\varphi}_j:\mathcal X\times[0,1)^{D-1}\to\mathcal X$ and the ordinary chain has independent $\bm{\Theta}_j\sim\operatorname{Unif}([0,1)^{D-1})$.

Suppose first that a scalar sorting key $\mathfrak h:\mathcal X\to\mathbb R$ has been chosen. At step $j$, order the old particles by this key, writing the resulting permutation as $\varsigma_j$, so that $\mathfrak h(\bm{X}_{\varsigma_j(0),j-1})\le\cdots\le \mathfrak h(\bm{X}_{\varsigma_j(N-1),j-1})$. Generate a freshly randomized point set $P_{N,j}=\{(V_{k,j},\bm{\Theta}_{k,j})\}_{k=0}^{N-1}\subset[0,1)^D$. The first coordinate $V_{k,j}$ serves only for sorting and matching; the remaining coordinates $\bm{\Theta}_{k,j}\in[0,1)^{D-1}$ drive the transition. Sort the complete points by their first coordinate, $V_{(0),j}\le\cdots\le V_{(N-1),j}$, retaining each point’s associated transition inputs, and match the particles and points by rank:

Fresh randomization is used at every transition, independently across time, while within a transition the low-discrepancy structure couples ordered particle ranks with transition inputs. The aim is to approximate the corresponding Markov-chain distribution more accurately than $N$ independently updated particles; in one dimension this can be viewed as improving the empirical CDF approximation.

The step function

records the old particle assigned to each rank interval. When the first coordinate is exactly $N$-stratified, $\bm{Z}_{N,j-1}(V_{(i),j})=\bm{X}_{\varsigma_j(i),j-1}$ almost surely, so evaluating the update through this step function gives precisely the rank-matched particle update. The scrambled-net condition used here is stated in ass:points.

The construction therefore requires an ordering of the particle states. For $p=1$, we use the natural numerical order, $\mathfrak h(x)=x$. For joint updates with $p>1$, there is no canonical natural ordering. Among the batch, split and Hilbert-based sorts reviewed in [12 Section 3], we use Hilbert ordering following [16]. Choose a componentwise continuous strictly increasing transformation $\bm{\psi}(\bm{x}) = (\psi_1(x_1),\ldots,\psi_p(x_p)): \mathcal X\to[0,1]^p,$ and write $\mathcal J:=\bm{\psi}(\mathcal X)\subseteq[0,1]^p.$ Thus $\bm{\psi}$ is a bijection from $\mathcal X$ onto $\mathcal J$, with componentwise inverse $\bm{\psi}^{-1}:\mathcal J\to\mathcal X$. Componentwise logistic or Gaussian-CDF transformations on unbounded state spaces are examples [16 Section 2.4]. Let $\bm{\mathscr H}_p:[0,1]\to[0,1]^p$ be a Hilbert curve and choose a Borel measurable right inverse $\iota_p:[0,1]^p\to[0,1]$ satisfying $\bm{\mathscr H}_p(\iota_p(\bm{z}))=\bm{z}$ [17 Proposition 2]. The scalar sorting key is $\mathfrak h(\bm{x})=\iota_p(\bm{\psi}(\bm{x}))$. In computation we use a finite-resolution Hilbert index [12 Section 3.2]. The transformation $\bm{\psi}$ is used only for sorting; the target, proposal and MH update remain defined on the original state space.

Array-RQMC MH and fixed-horizon results

The Array-RQMC-MH algorithm