Kernel Adaptive Metropolis-Hastings

Introduction

The choice of the proposal distribution is known to be crucial for the design of Metropolis-Hastings algorithms, and methods for adapting the proposal distribution to increase the sampler’s efficiency based on the history of the Markov chain have been widely studied. These methods often aim to learn the covariance structure of the target distribution, and adapt the proposal accordingly. Adaptive MCMC samplers were first studied by Haario et al. (1999, 2001), where the authors propose to update the proposal distribution along the sampling process. Based on the chain history, they estimate the covariance of the target distribution and construct a Gaussian proposal centered at the current chain state, with a particular choice of the scaling factor from Gelman et al. (1996). More sophisticated schemes are presented by Andrieu & Thoms (2008), e.g., adaptive scaling, component-wise scaling, and principal component updates.

While these strategies are beneficial for distributions that show high anisotropy (e.g., by ensuring the proposal uses the right scaling in all principal directions), they may still suffer from low acceptance probability and slow mixing when the target distributions are strongly nonlinear, and the directions of large variance depend on the current location of the sampler in the support. In the present work, we develop an adaptive Metropolis-Hastings algorithm in which samples are mapped to a reproducing kernel Hilbert space, and the proposal distribution is chosen according to the covariance in this feature space (Schölkopf et al., 1998; Smola et al., 2001). Unlike earlier adaptive approaches, the resulting proposal distributions are locally adaptive in input space, and oriented towards nearby regions of high density, rather than simply matching the global covariance structure of the distribution. Our approach combines a move in the feature space with a stochastic step towards the nearest input space point, where the feature space move can be analytically integrated out. Thus, the implementation of the procedure is straightforward: the proposal is simply a multivariate Gaussian in the input space, with location-dependent covariance which is informed by the feature space representation of the target. Furthermore, the resulting Metropolis-Hastings sampler only requires the ability to evaluate the unnormalized density of the target (or its unbiased estimate, as in Pseudo-Marginal MCMC of Andrieu & Roberts, 2009), and no gradient evaluation is needed, making it applicable to situations where more sophisticated schemes based on Hamiltonian Monte Carlo (HMC) or Metropolis Adjusted Langevin Algorithms (MALA) (Roberts & Stramer, 2003; Girolami & Calderhead, 2011) cannot be applied.

We begin our presentation in Section 2, with a brief overview of existing adaptive Metropolis approaches; we also review covariance operators in the RKHS. Based on these operators, we describe a sampling strategy for Gaussian measures in the RKHS in Section 3, and introduce a cost function for constructing proposal distributions. In Section 4, we outline our main algorithm, termed Kernel Adaptive Metropolis-Hastings (MCMC Kameleon). We provide experimental comparisons with other fixed and adaptive samplers in Section 5, where we show superior performance in the context of Pseudo-Marginal MCMC for Bayesian classification, and on synthetic target distributions with highly nonlinear shape.

Background

where $\nu=2.38/\sqrt{d}$ is a fixed scaling factor from Gelman et al. (1996). This choice of scaling factor was shown to be optimal (in terms of efficiency measures) for the usual Metropolis algorithm. While this optimality result does not hold for Adaptive Metropolis, it can nevertheless be used as a heuristic. Alternatively, the scale $\nu$ can also be adapted at each step as in Andrieu & Thoms (2008, Algorithm 4) to obtain the acceptance rate from Gelman et al. (1996), $a^{*}=0.234$.

RKHS Embeddings and Covariance Operators.

Our approach is based on the idea that the nonlinear support of a target density may be learned using Kernel Principal Component Analysis (Kernel PCA) (Schölkopf et al., 1998; Smola et al., 2001), this being linear PCA on the empirical covariance operator in the RKHS, $C_{\mathbf{z}}=\frac{1}{n}\sum_{i=1}^{n}k(\cdot,z_{i})\otimes k(\cdot,z_{i})-\mu_{\mathbf{z}}\otimes\mu_{\mathbf{z}}$, computed on the sample $\mathbf{z}$ defined above. The empirical covariance operator behaves as expected: applying the tensor product definition gives $\left\langle f,C_{\mathbf{z}}g\right\rangle_{\mathcal{H}_{k}}=\frac{1}{n}\sum_{i=1}^{n}f(z_{i})g(z_{i})-\left(\frac{1}{n}\sum_{i=1}^{n}f(z_{i})\right)\left(\frac{1}{n}\sum_{i=1}^{n}g(z_{i})\right)$. By analogy with algorithms which use linear PCA directions to inform M-H proposals (Andrieu & Thoms, 2008, Algorithm 8), nonlinear PCA directions can be encoded in the proposal construction, as described in Appendix C. Alternatively, one can focus on a Gaussian measure on the RKHS determined by the empirical covariance operator $C_{\mathbf{z}}$ rather than extracting its eigendirections, which is the approach we pursue in this contribution. This generalizes the proposal (1), which considers the Gaussian measure induced by the empirical covariance matrix on the original space.

Sampling in RKHS

We next describe the proposal distribution at iteration $t$ of the MCMC chain. We will assume that a subset of the chain history, denoted $\mathbf{z}=\left\{z_{i}\right\}_{i=1}^{n}$, $n\leq t-1$, is available. Our proposal is constructed by first considering the samples in the RKHS associated to the empirical covariance operator, and then performing a gradient descent step on a cost function associated with those samples.

We will work with the Gaussian measure on the RKHS $\mathcal{H}_{k}$ with mean $k(\cdot,y)$ and covariance $\nu^{2}C_{\mathbf{z}}$, where $\mathbf{z}=\left\{z_{i}\right\}_{i=1}^{n}$ is the subset of the chain history. While there is no analogue of a Lebesgue measure in an infinite dimensional RKHS, it is instructive (albeit with some abuse of notation) to denote this measure in the “density form” $\mathcal{N}(f\,;\,k(\cdot,y),\nu^{2}C_{\mathbf{z}})\propto\exp\left(-\frac{1}{2\nu^{2}}\left\langle f-k(\cdot,y),C_{\mathbf{z}}^{-1}(f-k(\cdot,y))\right\rangle_{\mathcal{H}_{k}}\right)$. As $C{}_{\mathbf{z}}$ is a finite-rank operator, this measure is supported only on a finite-dimensional affine space $k(\cdot,y)+\mathcal{H}_{\mathbf{z}}$, where $\mathcal{H}_{\mathbf{z}}=\textrm{span}\left\{k(\cdot,z_{i})\right\}_{i=1}^{n}$ is the subspace spanned by the canonical features of $\mathbf{z}$. It can be shown that a sample from this measure has the form $f=k(\cdot,y)+\sum_{i=1}^{n}\beta_{i}\left[k(\cdot,z_{i})-\mu_{\mathbf{z}}\right],$ where $\beta\sim\mathcal{N}(0,\frac{\nu^{2}}{n}I)$ is isotropic. Indeed, to see that $f$ has the correct covariance structure, note that:

Due to the equivalence in the RKHS between a Gaussian measure and a Gaussian Process (GP) (Berlinet & Thomas-Agnan, 2004, Ch. 4), we can think of the RKHS samples $f$ as trajectories of the GP with mean $m(x)=k(x,y)$ and covariance function

The covariance function $\kappa$ of this GP is therefore the kernel $k$ convolved with itself with respect to the empirical measure associated to the samples $\mathbf{z}$, and draws from this GP therefore lie in a smaller RKHS; see Saitoh (1997, p. 21) for details.

Obtaining Target Samples through Gradient Descent.

In general, this is a non-convex minimization problem, and may be difficult to solve (Bakir et al., 2003). Rather than solving it for every new vector of coefficients $\beta$, which would lead to an excessive computational burden for every proposal made, we simply make a single descent step along the gradient of the cost function,

where $\eta$ is a gradient step size parameter and $\xi\sim\mathcal{N}(0,\gamma^{2}I)$ is an additional isotropic ’exploration’ term after the gradient step. It will be useful to split the scaled gradient at $y$ into two terms as $\eta\nabla_{x}g(x)|_{x=y}=\eta\left(a_{y}-M_{\mathbf{z},y}H\beta\right)$, where $a_{y}=\nabla_{x}k(x,x)|_{x=y}-2\nabla_{x}k(x,y)|_{x=y}$,

is a $d\times n$ matrix, and $H=I-\frac{1}{n}\mathbf{1}_{n\times n}$ is the $n\times n$ centering matrix.

Figure 1 plots $g(x)$ and its gradients for several samples of $\beta$-coefficients, in the case where the underlying $\mathbf{z}$-samples are from the two-dimensional nonlinear Banana target distribution of Haario et al. (1999). It can be seen that $g$ may have multiple local minima, and that it varies most along the high-density regions of the Banana distribution.

MCMC Kameleon Algorithm

We now have a recipe to construct a proposal that is able to adapt to the local covariance structure for the current chain state $y$. This proposal depends on a subset of the chain history $\mathbf{z}$, and is denoted by $q_{\mathbf{z}}(\cdot|y)$. While we will later simplify this proposal by integrating out the moves in the RKHS, it is instructive to think of the proposal generating process as:

Sample $\beta\sim\mathcal{N}(0,\nu^{2}I)$ ($n\times 1$ normal of RKHS coefficients).

This represents an RKHS sample $f=k(\cdot,y)+\sum_{i=1}^{n}\beta_{i}\left[k(\cdot,z_{i})-\mu_{\mathbf{z}}\right]$ which is the goal of the cost function $g(x)$.

Move along the gradient of $g$: $x^{*}=y-\eta\nabla_{x}g(x)|_{x=y}+\xi.$

This gives a proposal $x^{*}|y,\beta\sim\mathcal{N}(y-\eta a_{y}+\eta M_{\mathbf{z},y}H\beta,\gamma^{2}I)$ ($d\times 1$ normal in the original space).

Our first step in the derivation of the explicit proposal density is to show that as long as $k$ is a differentiable positive definite kernel, the term $a_{y}$ vanishes.

Let $k$ be a differentiable positive definite kernel. Then $a_{y}=\nabla_{x}k(x,x)|_{x=y}-2\nabla_{x}k(x,y)|_{x=y}=0$.

Since $a_{y}=0$, the gradient step size $\eta$ always appears together with $\beta$, so we merge $\eta$ and the scale $\nu$ of the $\beta$-coefficients into a single scale parameter, and set $\eta=1$ henceforth. Furthermore, since both $p(\beta)$ and $p_{\mathbf{z}}(x^{*}|y,\beta)$ are multivariate Gaussian densities, the proposal density $q_{\mathbf{z}}(x^{*}|y)=\int p(\beta)p_{\mathbf{z}}(x^{*}|y,\beta)d\beta$ can be computed analytically. We therefore get the following closed form expression for the proposal distribution.

$q_{\mathbf{z}}(\cdot|y)=\mathcal{N}(y,\gamma^{2}I+\nu^{2}M_{\mathbf{z},y}HM_{\mathbf{z},y}^{\top})$.

Proofs of the above Propositions are given in Appendix A.

With the derived proposal distribution, we proceed with the standard Metropolis-Hastings accept/reject scheme, where the proposed sample $x^{*}$ is accepted with probability

giving rise to the MCMC Kameleon Algorithm. Note that each $\pi(x^{*})$ and $\pi(x_{t})$ could be replaced by their unbiased estimates without impacting the invariant distribution (Andrieu & Roberts, 2009).

The constructed family of proposals encodes local structure of the target distribution, which is learned based on the subsample $\mathbf{z}$. Figure 2 depicts the regions that contain 95% of the mass of the proposal distribution $q_{\mathbf{z}}(\cdot|y)$ at various states $y$ for a fixed subsample $\mathbf{z}$, where the Banana target is used (details in Section 5). More examples of proposal contours can be found in Appendix B.

2 Properties of the Algorithm

MCMC Kameleon requires a subsample $\mathbf{z}=\left\{z_{i}\right\}_{i=1}^{n}$ at each iteration of the algorithm, and the proposal distribution $q_{\mathbf{z}}(\cdot|y)$ is updated each time a new subsample $\mathbf{z}$ is obtained. It is well known that a chain which keeps adapting the proposal distribution need not converge to the correct target (Andrieu & Thoms, 2008). To guarantee convergence, we introduce adaptation probabilities $\left\{p_{t}\right\}_{t=0}^{\infty}$, such that $p_{t}\to 0$ and $\sum_{t=1}^{\infty}p_{t}=\infty$, and at iteration $t$ we update the subsample $\mathbf{z}$ with probability $p_{t}$. As adaptations occur with decreasing probability, Theorem 1 of Roberts & Rosenthal (2007) implies that the resulting algorithm is ergodic and converges to the correct target. Another straightforward way to guarantee convergence is to fix the set $\mathbf{z}=\left\{z_{i}\right\}_{i=1}^{n}$ after a “burn-in” phase; i.e., to stop adapting Roberts & Rosenthal (2007, Proposition 2). In this case, a “burn-in” phase is used to get a rough sketch of the shape of the distribution: the initial samples need not come from a converged or even valid MCMC chain, and it suffices to have a scheme with good exploratory properties, e.g., Welling & Teh (2011). In MCMC Kameleon, the term $\gamma$ allows exploration in the initial iterations of the chain (while the subsample $\mathbf{z}$ is still not informative about the structure of the target) and provides regularization of the proposal covariance in cases where it might become ill-conditioned. Intuitively, a good approach to setting $\gamma$ is to slowly decrease it with each adaptation, such that the learned covariance progressively dominates the proposal.

Symmetry of the proposal.

In Haario et al. (2001), the proposal distribution is asymptotically symmetric due to the vanishing adaptation property. Therefore, the authors compute the standard Metropolis acceptance probability. In our case, the proposal distribution is a Gaussian with mean at the current state of the chain $x_{t}=y$ and covariance $\gamma^{2}I+\nu^{2}M_{\mathbf{z},y}HM_{\mathbf{z},y}^{\top}$, where $M_{\mathbf{z},y}$ depends both on the current state $y$ and a random subsample $\mathbf{z}=\left\{z_{i}\right\}_{i=1}^{n}$ of the chain history $\left\{x_{i}\right\}_{i=0}^{t-1}$. This proposal distribution is never symmetric (as covariance of the proposal always depends on the current state of the chain), and therefore we use the Metropolis-Hastings acceptance probability to reflect this.

Relationship to MALA and Manifold MALA.

The Metropolis Adjusted Langevin Algorithm (MALA) algorithm uses information about the gradient of the log-target density at the current chain state to construct a proposed point for the Metropolis step. Our approach does not require that the log-target density gradient be available or computable. Kernel gradients in the matrix $M_{\mathbf{z},y}$ are easily obtained for commonly used kernels, including the Gaussian kernel (see section 4.3), for which the computational complexity is equal to evaluating the kernel itself. Moreover, while standard MALA simply shifts the mean of the proposal distribution along the gradient and then adds an isotropic exploration term, our proposal is centered at the current state, and it is the covariance structure of the proposal distribution that coerces the proposed points to belong to the high-density regions of the target. It would be straightforward to modify our approach to include a drift term along the gradient of the log-density, should such information be available, but it is unclear whether this would provide additional performance gains. Further work is required to elucidate possible connections between our approach and the use of a preconditioning matrix (Roberts & Stramer, 2003) in the MALA proposal; i.e., where the exploration term is scaled with appropriate metric tensor information, as in Riemannian manifold MALA (Girolami & Calderhead, 2011).

3 Examples of Covariance Structure for Standard Kernels

The proposal distributions in MCMC Kameleon are dependant on the choice of the kernel $k$. To gain intuition regarding their covariance structure, we give two examples below.

In the case of a linear kernel $k(x,x^{\prime})=x^{\top}x^{\prime}$, we obtain $M_{\mathbf{z},y}=2\left[\nabla_{x}x^{\top}z_{1}|_{x=y},\ldots,\nabla_{x}x^{\top}z_{n}|_{x=y}\right]=2\mathbf{Z}^{\top}$, so the proposal is given by $q_{\mathbf{z}}(\cdot|y)=\mathcal{N}(y,\gamma^{2}I+4\nu^{2}\mathbf{Z}^{\top}H\mathbf{Z})$; thus, the proposal simply uses the scaled empirical covariance $\mathbf{Z}^{\top}H\mathbf{Z}$ just like standard Adaptive Metropolis (Haario et al., 1999), with an additional isotropic exploration component, and depends on $y$ only through the mean.

Gaussian kernel.

In the case of a Gaussian kernel $k(x,x^{\prime})=\exp\left(-\frac{\left\|x-x^{\prime}\right\|_{2}^{2}}{2\sigma^{2}}\right)$, since $\nabla_{x}k(x,x^{\prime})=\frac{1}{\sigma^{2}}k(x,x^{\prime})(x^{\prime}-x)$, we obtain

Consider how this encodes the covariance structure of the target distribution:

As the first two terms dominate, the previous points $z_{a}$ which are close to the current state $y$ (for which $k(y,z_{a})$ is large) have larger weights, and thus they have more influence in determining the covariance of the proposal at $y$.

Matérn kernel.

In the Matérn family of kernels $k_{\vartheta,\rho}(x,x^{\prime})=\frac{2^{1-\vartheta}}{\Gamma(\vartheta)}\left(\frac{\left\|x-x^{\prime}\right\|_{2}}{\rho}\right)^{\vartheta}K_{\vartheta}\left(\frac{\left\|x-x^{\prime}\right\|_{2}}{\rho}\right)$, where $K_{\vartheta}$ is the modified Bessel function of the second kind, we obtain a form of the covariance structure very similar to that of the Gaussain kernel. In this case, $\nabla_{x}k_{\vartheta,\rho}(x,x^{\prime})=\frac{1}{2\rho^{2}(\vartheta-1)}k_{\vartheta-1,\rho}(x,x^{\prime})(x^{\prime}-x)$, so the only difference (apart from the scalings) to (5) is that the weights are now determined by a “rougher” kernel $k_{\vartheta-1,\rho}$ of the same family.

Experiments

In the experiments, we compare the following samplers: (SM) Standard Metropolis with the isotropic proposal $q(\cdot|y)=\mathcal{N}(y,\nu^{2}I)$ and scaling $\nu=2.38/\sqrt{d}$, (AM-FS) Adaptive Metropolis with a learned covariance matrix and fixed scaling $\nu=2.38/\sqrt{d}$, (AM-LS) Adaptive Metropolis with a learned covariance matrix and scaling learned to bring the acceptance rate close to $\alpha^{*}=0.234$ as described in Andrieu & Thoms (2008, Algorithm 4), and (KAMH-LS) MCMC Kameleon with the scaling $\nu$ learned in the same fashion ($\gamma$ was fixed to 0.2), and which also stops adapting the proposal after the burn-in of the chain (in all experiments, we use a random subsample $\mathbf{z}$ of size $n=1000$, and a Gaussian kernel with bandwidth selected according to the median heuristic). We consider the following nonlinear targets: (1) the posterior distribution of Gaussian Process (GP) classification hyperparameters (Filippone & Girolami, 2014) on the UCI glass dataset, and (2) the synthetic banana-shaped distribution of Haario et al. (1999) and a flower-shaped disribution concentrated on a circle with a periodic perturbation.

In the first experiment, we illustrate usefulness of the MCMC Kameleon sampler in the context of Bayesian classification with GPs (Williams & Barber, 1998). Consider the joint distribution of latent variables $\mathbf{f}$, labels $\mathbf{y}$ (with covariate matrix $X$), and hyperparameters $\theta$, given by

where $\left\{\mathbf{f}^{(i)}\right\}_{i=1}^{n_{\textrm{imp}}}\sim q(\mathbf{f}|\theta)$ are $n_{imp}$ importance samples. In Filippone & Girolami (2014), the importance distribution $q(\mathbf{f}|\theta)$ is chosen as the Laplacian or as the Expectation Propagation (EP) approximation of $p(\mathbf{f}|\mathbf{y},\theta)\propto p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\theta)$, leading to state-of-the-art results.

We consider the UCI Glass dataset (Bache & Lichman, 2013), where classification of window against non-window glass is sought. Due to the heterogeneous structure of each of the classes (i.e., non-window glass consists of containers, tableware and headlamps), there is no single consistent set of lengthscales determining the decision boundary, so one expects the posterior of the covariance bandwidths $\theta_{d}$ to have a complicated (nonlinear) shape. This is illustrated by the plot of the posterior projections to the dimensions 2 and 7 (out of 9) in Figure 3. Since the ground truth for the hyperparameter posterior is not available, we initially ran 30 Standard Metropolis chains for 500,000 iterations (with a 100,000 burn-in), kept every 1000-th sample in each of the chains, and combined them. The resulting samples were used as a benchmark, to evaluate the performance of shorter single-chain runs of SM, AM-FS, AM-LS and KAMH-LS. Each of these algorithms was run for 100,000 iterations (with a 20,000 burnin) and every 20-th sample was kept. Two metrics were used in evaluating the performance of the four samplers, relative to the large-scale benchmark. First, the distance $\left\|\hat{\mu}_{\theta}-\mu_{\theta}^{b}\right\|_{2}$ was computed between the mean $\hat{\mu}_{\theta}$ estimated from each of the four sampler outputs, and the mean $\mu_{\theta}^{b}$ on the benchmark sample (Fig. 4, left), as a function of sample size. Second, the MMD (Borgwardt et al., 2006; Gretton et al., 2007) was computed between each sampler output and the benchmark sample, using the polynomial kernel $\left(1+\left\langle\theta,\theta^{\prime}\right\rangle\right)^{3}$; i.e., the comparison was made in terms of all mixed moments of order up to 3 (Fig. 4, right). The figures indicate that KAMH-LS approximates the benchmark sample better than the competing approaches, where the effect is especially pronounced in the high order moments, indicating that KAMH-LS thoroughly explores the distribution support in a relatively small number of samples.

We emphasise that, as for any pseudo-marginal MCMC scheme, neither the likelihood itself, nor any higher-order information about the marginal posterior target $p(\theta|\mathbf{y})$, are available. This makes HMC or MALA based approaches such as (Roberts & Stramer, 2003; Girolami & Calderhead, 2011) unsuitable for this problem, so it is very difficult to deal with strongly nonlinear posterior targets. In contrast, as indicated in this example, the MCMC Kameleon scheme is able to effectively sample from such nonlinear targets, and outperforms the vanilla Metropolis methods, which are the only competing choices in the pseudo-marginal context.

In addition, since the bulk of the cost for pseudo-marginal MCMC is in importance sampling in order to obtain the acceptance ratio, the additional cost imposed by KAMH-LS is negligible. Indeed, we observed that there is an increase of only 2-3% in terms of effective computation time in comparison to all other samplers, for the chosen size of the chain history subsample ($n=1000$).

2 Synthetic examples

Flower target.

The second target distribution we consider is the $d$-dimensional flower target $\mathcal{F}(r_{0},A,\omega,\sigma)$, with

This distribution concentrates around the $r_{0}$-circle with a periodic perturbation (with amplitude $A$ and frequency $\omega$) in the first two dimensions.

In these examples, exact quantile regions of the targets can be computed analytically, so we can directly assess performance without the need to estimate distribution distances on the basis of samples (i.e., by estimating MMD to the benchmark sample). We compute the following measures of performance (similarly as in Haario et al. (1999); Andrieu & Thoms (2008)) based on the chain after burn-in: average acceptance rate, norm of the empirical mean (the true mean is by construction zero for all targets), and the deviation of the empirical quantiles from the true quantiles. We consider 8-dimensional target distributions: the moderately twisted $\mathcal{B}(0.03,100)$ banana target (Figure 5, top) and the strongly twisted $\mathcal{B}(0.1,100)$ banana target (Figure 5, middle) and $\mathcal{F}(10,6,6,1)$ flower target (Figure 5, bottom).

The results show that MCMC Kameleon is superior to the competing samplers. Since the covariance of the proposal adapts to the local structure of the target at the current chain state, as illustrated in Figure 2, MCMC Kameleon does not suffer from wrongly scaled proposal distributions. The result is a significantly improved quantile performance in comparison to all competing samplers, as well as a comparable or superior norm of the empirical mean. SM has a significantly larger norm of the empirical mean, due to its purely random walk behavior (e.g., the chain tends to get stuck in one part of the space, and is not able to traverse both tails of the banana target equally well). AM with fixed scale has a low acceptance rate (indicating that the scaling of the proposal is too large), and even though the norm of the empirical mean is much closer to the true value, quantile performance of the chain is poor. Even if the estimated covariance matrix closely resembles the true global covariance matrix of the target, using it to construct proposal distributions at every state of the chain may not be the best choice. For example, AM correctly captures scalings along individual dimensions for the flower target (the norm of its empirical mean is close to its true value of zero) but fails to capture local dependence structure. The flower target, due to its symmetry, has an isotropic covariance in the first two dimensions – even though they are highly dependent. This leads to a mismatch in the scale of the covariance and the scale of the target, which concentrates on a thin band in the joint space. AM-LS has the “correct” acceptance rate, but the quantile performance is even worse, as the scaling now becomes too small to traverse high-density regions of the target.

Conclusions

We have constructed a simple, versatile, adaptive, gradient-free MCMC sampler that constructs a family of proposal distributions based on the sample history of the chain. These proposal distributions automatically conform to the local covariance structure of the target distribution at the current chain state. In experiments, the sampler outperforms existing approaches on nonlinear target distributions, both by exploring the entire support of these distributions, and by returning accurate empirical quantiles, indicating faster mixing. Possible extensions include incorporating additional parametric information about the target densities, and exploring the tradeoff between the degree of sub-sampling of the chain history and convergence of the sampler.

Python implementation of MCMC Kameleon is available at https://github.com/karlnapf/kameleon-mcmc.

Acknowledgments.

D.S., H.S., M.L.G. and A.G. acknowledge support of the Gatsby Charitable Foundation. We thank Mark Girolami for insightful discussions and the anonymous reviewers for useful comments.

References

Appendix A Proofs

Let $k$ be a differentiable positive definite kernel. Then $\nabla_{x}k(x,x)|_{x=y}-2\nabla_{x}k(x,y)|_{x=y}=0$.

Next, we consider the map $\kappa_{y}(x)=k(x,y)=\left\langle\varphi(x),\varphi(y)\right\rangle_{\mathcal{H}}$, i.e., $\kappa_{y}=\psi_{y}\circ\varphi$ where $\psi_{y}(f)=\left\langle f,\varphi(y)\right\rangle_{\mathcal{H}}$. Since $\psi_{y}$ is a linear scalar function on $\mathcal{H}$, $D\psi_{y}\left(f\right)=\left\langle\varphi(y),\cdot\right\rangle_{\mathcal{H}}$. Again, by the chain rule:

$q_{\mathbf{z}}(\cdot|y)=\mathcal{N}(y,\gamma^{2}I+\nu^{2}M_{\mathbf{z},y}HM_{\mathbf{z},y}^{\top})$.

and application of the standard Gaussian integral

This is just a $d$-dimensional Gaussian density where both the mean and covariance will, in general, depend on $y$. Let us consider the exponent

where $R^{-1}=\frac{1}{\gamma^{2}}(I-\frac{1}{\gamma^{2}}M_{\mathbf{z},y}H\Sigma HM_{\mathbf{z},y}^{\top})$. We can simplify the covariance $R$ using the Woodbury identity to obtain:

Therefore, the proposal density is $q_{\mathbf{z}}(\cdot|y)=\mathcal{N}(y,\gamma^{2}I+\nu^{2}M_{\mathbf{z},y}HM_{\mathbf{z},y}^{\top})$. ∎

Appendix B Further details on synthetic experiments

The $d$-dimensional flower target $\mathcal{F}(r_{0},A,\omega,\sigma)$ is given by

This distribution concentrates around the $r_{0}$-circle with a periodic perturbation (with amplitude $A$ and frequency $\omega$) in the first two dimensions. For $A=0$, we obtain a band around the $r_{0}$-circle, which we term the ring target. Figure 6 gives the contour plots of the MCMC Kameleon proposal distributions on two instances of the flower target.

Convergence statistics for the Banana target.

Figure 7 illustrates how the norm of the mean and quantile deviation (shown for $0.5$-quantile) for the strongly twisted Banana target decrease as a function of the number of iterations. This shows that the trends observed in the main text persist along the evolution of the whole chain.

Appendix C Principal Components Proposals

An alternative approach to the standard adaptive Metropolis, discussed in Andrieu & Thoms (2008, Algorithm 8), is to extract $m\leq d$ principal eigenvalue-eigenvector pairs $\left\{\left(\lambda_{j},v_{j}\right)\right\}_{j=1}^{m}$ from the estimated covariance matrix $\Sigma_{\mathbf{z}}$ and use the proposal that takes form of a mixture of one-dimensional random walks along the principal eigendirections

In other words, given the current chain state $y$, the $j$-th principal eigendirection is chosen with probability $\omega_{j}$ (choice $\omega_{j}=\lambda_{j}/\sum_{l=1}^{m}\lambda_{l}$ is suggested), and the proposed point is

with $\rho\sim\mathcal{N}(0,1)$. Note that each eigendirection may have a different scaling factor $\nu_{j}$ in addition to the scaling with the eigenvalue.

We can consider an analogous version of the update (8) performed in the RKHS

Now, we can construct the MCMC PCA-Kameleon by simply substituting $\beta$-coefficients with $\rho\nu_{j}\alpha^{(j)}$, where $j$ is the selected eigendirection, and $\nu_{j}$ is the scaling factor associated to the $j$-th eigendirection. We have the following steps:

Perform eigendecomposition of $HKH$ to obtain the $m\leq n$ eigenvectors $\left\{\alpha_{j}\right\}_{j=1}^{m}.$

Draw $j\sim\text{{Discrete}}\left[\omega_{1},\ldots,\omega_{m}\right]$

$x^{*}|y,\rho,j\sim\mathcal{N}(y+\rho\nu_{j}M_{\mathbf{z},y}H\alpha^{(j)},\gamma^{2}I)$ ($d\times 1$ normal in the original space)

Similarly as before, we can simplify the proposal by integrating out the scale $\rho$ of the moves in the RKHS.

$q_{\mathbf{z}}(\cdot|y)=\sum_{j=1}^{m}\omega_{j}\mathcal{N}(y,\gamma^{2}I+\nu_{j}^{2}M_{\mathbf{z},y}H\alpha^{(j)}\left(\alpha^{(j)}\right)^{\top}HM_{\mathbf{z},y}^{\top})$.

where $R^{-1}=\frac{1}{\gamma^{2}}\left(I-\frac{\nu_{j}^{2}\sigma^{2}}{\gamma^{2}}M_{\mathbf{z},y}H\alpha^{(j)}\left(\alpha^{(j)}\right)^{\top}HM_{\mathbf{z},y}^{\top}\right)$. We can simplify the covariance $R$ using the Woodbury identity to obtain:

The claim follows after summing over the choice $j$ of the eigendirection (w.p. $\omega_{j}$).∎