Can the adaptive Metropolis algorithm collapse without the covariance lower bound?
Full text
This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Can the adaptive Metropolis algorithm collapse without the covariance lower bound? Vihola, Matti Vihola, M. (2011). Can the adaptive Metropolis algorithm collapse without the covariance lower bound?. Electronic Journal of Probability, 16, 45-75. DOI:10.1214/EJP.v16-840. 2011
Electronic Journal of Probability Vol.16 (2011), Paper no. 2, pages 45–75. Journal URL http://www.math.washington.edu/~ejpecp/ Can the adaptive Metropolis algorithm collapse without the covariance lower bound?∗ Matti Vihola Department of Mathematics and Statistics University of Jyväskylä P.O.Box 35 FI-40014 University of Jyväskyla Finland [email protected] http://iki.fi/mvihola/ Abstract The Adaptive Metropolis (AM) algorithm is based on the symmetric random-walk Metropolis algorithm. The proposal distribution has the following time-dependent covariance matrix at step n+1 Sn=Cov(X1,...,Xn) + εI, that is, the sample covariance matrix of the history of the chain plus a (small) constant ε > 0 multiple of the identity matrix I. The lower bound on the eigenvalues of Sninduced by the factor εIis theoretically convenient, but practically cumbersome, as a good value for the parameter εmay not always be easy to choose. This article considers variants of the AM algorithm that do not explicitly bound the eigenvalues of Snaway from zero. The behaviour of Snis studied in detail, indicating that the eigenvalues of Sndo not tend to collapse to zero in general. In dimension one, it is shown that Snis bounded away from zero if the logarithmic target density is uniformly continuous. For a modification of the AM algorithm including an additional fixed ∗The author was supported by the Academy of Finland, projects no. 110599 and 201392, by the Finnish Academy of Science and Letters, Vilho, Yrjö and Kalle Väisälä Foundation, by the Finnish Centre of Excellence in Analysis and Dynamics Research, and by the Finnish Graduate School in Stochastics and Statistics 45
component in the proposal distribution, the eigenvalues of Snare shown to stay away from zero with a practically non-restrictive condition. This result implies a strong law of large numbers for super-exponentially decaying target distributions with regular contours. Key words: Adaptive Markov chain Monte Carlo, Metropolis algorithm, stability, stochastic approximation. AMS 2000 Subject Classification: Primary 65C40. Submitted to EJP on January 14, 2010, final version accepted November 27, 2011. 46
1 Introduction Adaptive Markov chain Monte Carlo (MCMC) methods have attracted increasing interest in the last few years, after the original work of Haario, Saksman, and Tamminen Haario et al. (2001) and the subsequent advances in the field Andrieu and Moulines (2006); Andrieu and Robert (2001); Atchadé and Rosenthal (2005); Roberts and Rosenthal (2007); see also the recent review Andrieu and Thoms (2008). Several adaptive MCMC algorithms have been proposed up to date, but the seminal Adaptive Metropolis (AM) algorithm Haario et al. (2001) is still one of the most applied methods, perhaps due to its simplicity and generality. The AM algorithm is a symmetric random-walk Metropolis algorithm, with an adaptive proposal distribution. The algorithm starts1at some point X1≡x1∈Rdwith an initial positive definite covariance matrix S1≡s1∈Rd×dand follows the recursion (S1) Let Yn+1=Xn+θS1/2 nWn+1, where Wn+1is an independent standard Gaussian random vector and θ > 0 is a constant. (S2) Accept Yn+1with probability min1, π(Yn+1) π(Xn)and let Xn+1=Yn+1; otherwise reject Yn+1and let Xn+1=Xn. (S3) Set Sn+1= Γ(X1,...,Xn+1). In the original work Haario et al. (2001) the covariance parameter is computed by Γ(X1,...,Xn+1) = 1 n n+1 X k=1 (Xk−Xn+1)(Xk−Xn+1)T+εI, (1) where Xn:=n−1Pn k=1Xkstands for the mean. That is, Sn+1is a covariance estimate of the history of the ‘Metropolis chain’ X1,...,Xn+1plus a small ε > 0 multiple of the identity matrix I∈Rd×d. The authors prove a strong law of large numbers (SLLN) for the algorithm, that is, n−1Pn k=1f(Xk)→RRdf(x)π(x)dxalmost surely as n→ ∞ for any bounded functional fwhen the target distribution πis bounded and compactly supported. Recently, SLLN was shown to hold also for πwith unbounded support, having super-exponentially decaying tails with regular contours and fgrowing at most exponentially in the tails Saksman and Vihola (2010). This article considers the original AM algorithm (S1)–(S3), without the lower bound induced by the factor εI. The proposal covariance function Γ, defined precisely in Section 2, is a consistent covariance estimator first proposed in Andrieu and Robert (2001). A special case of this estimator behaves asymptotically like the sample covariance in (1). Previous results indicate that if this algorithm is modified by truncating the eigenvalues of Snwithin explicit lower and upper bounds, the algorithm can be verified in a fairly general setting Atchadé and Fort (2010); Roberts and Rosenthal (2007). It is also possible to determine an increasing sequence of truncation sets for Sn, and modify the algorithm to include a re-projection scheme in order to verify the validity of the algorithm Andrieu and Moulines (2006). While technically convenient, such pre-defined bounds on the adapted covariance matrix Sncan be inconvenient in practice. Ill-defined values can affect the efficiency of the adaptive scheme dramatically, rendering the algorithm useless in the worst case. In particular, if the factor ε > 0 in the AM 1The initial ‘burn-in’ phase included in the original algorithm is not considered here. 47
algorithm is selected too large, the smallest eigenvalue of the true covariance matrix of πmay be well smaller than ε > 0, and the chain Xnis likely to mix poorly. Even though the re-projection scheme of Andrieu and Moulines (2006) avoids such behaviour by increasing truncation sets, which eventually contain the desirable values of the adaptation parameter, the practical efficiency of the algorithm is still strongly affected by the choice of these sets Andrieu and Thoms (2008). Without a lower bound on the eigenvalues of Sn(or a re-projection scheme), there is a potential danger of the covariance parameter Sncollapsing to singularity. In such a case, the increments Xn−Xn−1would be smaller and smaller, and the Xnchain could eventually get ‘stuck’. The empirical evidence suggests that this does not tend to happen in practice. The present results validate the empirical findings by excluding such a behaviour in different settings. After defining precisely the algorithms in Section 2, the above mentioned unconstrained AM algorithm is analysed in Section 3. First, the AM algorithm run on an improper uniform target π≡c>0 is studied. In such a case, the asymptotic expected growth rate of Snis characterised quite precisely, being e2θpnfor the original AM algorithm Haario et al. (2001). The behaviour of the AM algorithm in the uniform target setting is believed to be similar as in a situation where Snis small and the target πis smooth whence locally constant. The results support the strategy of choosing a ‘small’ initial covariance s1in practice, and letting the adaptation take care of expanding it to the proper size. In Section 3, it is also shown that in a one-dimensional setting and with a uniformly continuous logπ, the variance parameter Snis bounded away from zero. This fact is shown to imply, with the results in Saksman and Vihola (2010), a SLLN in the particular case of a Laplace target distribution. While this result has little practical value in its own right, it is the first case where the unconstrained AM algorithm is shown to preserve the correct ergodic properties. It shows that the algorithm possesses self-stabilising properties and further strengthens the belief that the algorithm would be stable and ergodic under a more general setting. Section 4 considers a slightly different variant of the AM algorithm, due to Roberts and Rosenthal Roberts and Rosenthal (2009), replacing (S1) with (S1’) With probability β, let Yn+1=Xn+Vn+1where Vn+1is an independent sample of qfix; otherwise, let Yn+1=Xn+θS1/2 nWn+1as in (S1). While omitting the parameter ε > 0, the proposal strategy (S1’) includes two additional parameters: the mixing probability β∈(0,1)and the fixed symmetric proposal distribution qfix. It has the advantage that the ‘worst case scenario’ having ill-defined qfix only ‘wastes’ the fixed proportion βof samples, while Sncan take any positive definite value on adaptation. This approach is analysed also in the recent preprint Bai et al. (2008), relying on a technical assumption that ultimately implies that Xnis bounded in probability. In particular, the authors show that if qfix is a uniform density on a ball having a large enough radius, then the algorithm is ergodic. Section 4 uses a perhaps more transparent argument to show that the proposal strategy (S1’) with a mild additional condition implies a sequence Snwith eigenvalues bounded away from zero. This fact implies a SLLN using the technique of Saksman and Vihola (2010), as shown in the end of Section 4. 48
2 The general algorithm Let us define a Markov chain (Xn,Mn,Sn)n≥1evolving in space Rd×Rd×Cdwith the state space Rdand Cd⊂Rd×dstanding for the positive definite matrices. The chain starts at an initial position X1≡x1∈Rd, with an initial mean2M1≡m1∈Rdand an initial covariance matrix S1≡s1∈ Cd. For n≥1, the chain is defined through the recursion Xn+1∼PqSn(Xn,·)(2) Mn+1:= (1−ηn+1)Mn+ηn+1Xn+1(3) Sn+1:= (1−ηn+1)Sn+ηn+1(Xn+1−Mn)(Xn+1−Mn)T. (4) Denoting the natural filtration of the chain as Fn:=σ(Xk,Mk,Sk: 1 ≤k≤n), the notation in (2) reads that PXn+1∈AFn=PqSn(Xn,A)for any measurable A⊂Rd. The Metropolis transition kernel Pqis defined for any symmetric probability density q(x,y) = q(x−y)through Pq(x,A):= 1 A(x)1−Zmin1, π(y) π(x)q(y−x)dy +ZA min1, π(y) π(x)q(y−x)dy where 1 Astands for the characteristic function of the set A. The proposal densities {qs}s∈Cdare defined as a mixture qs(z):= (1−β)˜ qs(z) + βqfix(z)(5) where the mixing constant β∈[0,1)determines the portion how often a fixed proposal density qfix is used instead of the adaptive proposal ˜ qs(z):=det(θs)−1/2˜ q(θ−1/2s−1/2z)with ˜ qbeing a ‘template’ probability density. Finally, the adaptation weights (ηn)n≥2⊂(0,1)appearing in (3) and (4) is assumed to decay to zero. One can verify that for β=0 this setting corresponds to the algorithm (S1)–(S3) of Section 1 with Wn+1having distribution ˜ q, and for β∈(0,1), (S1’) applies instead of (S1). Notice also that the original AM algorithm essentially fits this setting, with ηn:=n−1,β:=0 and if ˜ qsis defined slightly differently, being a Gaussian density with mean zero and covariance s+εI. Moreover, if one sets β=1, the above setting reduces to a non-adaptive symmetric random walk Metropolis algorithm with the increment proposal distribution qfix. 3 The unconstrained AM algorithm 3.1 Overview of the results This section deals with the unconstrained AM algorithm, that is, the algorithm described in Section 2 with the mixing constant β=0 in (5). Sections 3.2 and 3.3 consider the case of an improper uniform target distribution π≡cfor some constant c>0. This implies that (almost) every proposed sample is accepted and the recursion (2) reduces to Xn+1=Xn+θS1/2 nWn+1(6) 2A customary choice is to set m1=x1. 49
100102104106 −10 −5 0 5 ESn n Figure 1: An example of the exact development of ESn, when s1=1 and θ=0.01. The sequence (ESn)n≥1decreases until nis over 27,000 and exceeds the initial value only with nover 750,000. where (Wn)n≥2are independent realisations of the distribution ˜ q. Throughout this subsection, let us assume that the template proposal distribution ˜ qis spherically symmetric and the weight sequence is defined as ηn:=cn−γfor some constants c∈(0,1]and γ∈(1/2,1]. The first result characterises the expected behaviour of Snwhen (Xn)n≥2follows (6). Theorem 1. Suppose (Xn)n≥2follows the ‘adaptive random walk’ recursion (6), with EWnWT n=I. Then, for all λ > 1there is n0≥m such that for all n ≥n0and k ≥1, the following bounds hold 1 λ θ n+k X j=n+1pηj ≤logESn+k ESn≤λ θ n+k X j=n+1pηj . Proof. Theorem 1 is a special case of Theorem 12 in Section 3.2. Remark 2.Theorem 1 implies that with the choice ηn:=cn−γfor some c∈(0,1)and γ∈(1/2,1], the expectation grows with the speed ESn≃exp θpc 1−γ 2 n1−γ 2!. Remark 3.In the original setting (Haario et al. 2001) the weights are defined as ηn:=n−1and Theorem 1 implies that the asymptotic growth rate of E[Sn]is e2θpnwhen (Xn)n≥2follows (6). Suppose the value of Snis very small compared to the scale of a smooth target distribution π. Then, it is expected that most of the proposal are accepted, Xnbehaves almost as (6), and Snis expected to grow approximately at the rate e2θpnuntil it reaches the correct magnitude. On the other hand, simple deterministic bound implies that Sncan decay slowly, only with the polynomial speed n−1. Therefore, it may be safer to choose the initial s1small. Remark 4.The selection of the scaling parameter θ > 0 in the AM algorithm does not seem to affect the expected asymptotic behaviour Sndramatically. However, the choice 0 < θ ≪1 can result in an significant initial ‘dip’ of the adapted covariance values, as exemplified in Figure 1. Therefore, the values θ≪1 are to be used with care. In this case, the significance of a successful burn-in is also emphasised. 50
It may seem that Theorem 1 would automatically also ensure that Sn→ ∞ also path-wise. This is not, however, the case. For example, consider the probability space [0,1]with the Borel σ-algebra and the Lebesgue measure. Then (Mn,Fn)n≥1defined as Mn:=22n 1 [0,2−n)and Fn:=σ(Xk: 1 ≤ k≤n)is, in fact, a submartingale. Moreover, EMn=2n→∞, but Mn→0 almost surely. The AM process, however, does produce an unbounded sequence Sn. Theorem 5. Assume that (Xn)n≥2follows the ‘adaptive random walk’ recursion (6). Then, for any unit vector u ∈Rd, the process uTSnu→∞ almost surely. Proof. Theorem 5 is a special case of Theorem 18 in Section 3.3. In a one-dimensional setting, and when logπis uniformly continuous, the AM process can be approximated with the ‘adaptive random walk’ above, whenever Snis small enough. This yields Theorem 6. Assume d =1and logπis uniformly continuous. Then, there is a constant b >0such that liminfn→∞ Sn≥b. Proof. Theorem 6 is a special case of Theorem 18 in Section 3.4. Finally, having Theorem 6, it is possible to establish Theorem 7. Assume ˜ q is Gaussian, the one-dimensional target distribution is standard Laplace π(x):= 1 2e−|x|and the functional f :R→Rsatisfies supxe−γ|x||f(x)|<∞for some γ∈(0,1/2). Then, n−1Pn k=1f(Xk)→Rf(x)π(x)dx almost surely as n →∞. Proof. Theorem 7 is a special case of Theorem 21 in Section 3.4. Remark 8.In the case ηn:=n−1, Theorem 7 implies that the parameters Mnand Snof the adaptive chain converge to 0 and 2, that is, the true mean and variance of the target distribution π, respectively. Remark 9.Theorem 6 (and Theorem 7) could probably be extended to cover also targets πwith compact supports. Such an extension would, however, require specific handling of the boundary effects, which can lead to technicalities. 3.2 Uniform target: expected growth rate Define the following matrix quantities an:=E(Xn−Mn−1)(Xn−Mn−1)T(7) bn:=ESn(8) for n≥1, with the convention that a1≡0∈Rd×d. One may write using (3) and (6) Xn+1−Mn=Xn−Mn+θS1/2 nWn+1= (1−ηn)(Xn−Mn−1) + θS1/2 nWn+1. If EWnWT n=I, one may easily compute E(Xn+1−Mn)(Xn+1−Mn)T =1−ηn2E(Xn−Mn−1)(Xn−Mn−1)T+θ2ESn 51
since Wn+1is independent of Fnand zero-mean due to the symmetry of ˜ q. The values of (an)n≥2 and (bn)n≥2are therefore determined by the joint recursion an+1= (1−ηn)2an+θ2bn(9) bn+1= (1−ηn+1)bn+ηn+1an+1. (10) Observe that for any constant unit vector u∈Rd, the recursions (9) and (10) hold also for a(u) n+1:=EuT(Xn+1−Mn)(Xn+1−Mn)Tu b(u) n+1:=EuTSn+1u. The rest of this section therefore dedicates to the analysis if the one-dimensional recursions (9) and (10), that is, an,bn∈R+for all n≥1. The first result shows that the tail of (bn)n≥1is increasing. Lemma 10. Let n0≥1and suppose an0≥0, bn0>0and for n ≥n0the sequences anand bnfollow the recursions (9) and (10), respectively. Then, there is a m0≥n0such that (bn)n≥m0is strictly increasing. Proof. If θ≥1, we may estimate an+1≥(1−ηn)2an+bnimplying bn+1≥bn+ηn+1(1−ηn)2an for all n≥n0. Since bn>0 by construction, and therefore also an+1≥θ2bn>0, we have that bn+1>bnfor all n≥n0+1. Suppose then θ < 1. Solving an+1from (10) yields an+1=η−1 n+1bn+1−bn+bn Substituting this into (9), we obtain for n≥n0+1 η−1 n+1bn+1−bn+bn= (1−ηn)2η−1 nbn−bn−1+bn−1+θ2bn After some algebraic manipulation, this is equivalent to bn+1−bn=ηn+1 ηn (1−ηn)3(bn−bn−1) + ηn+1(1−ηn)2−1+θ2bn. (11) Now, since ηn→0, we have that (1−ηn)2−1+θ2>0 whenever nis greater than some n1. So, if we have for some n′>n1that bn′−bn′−1≥0, the sequence (bn)n≥n′is strictly increasing after n′. Suppose conversely that bn+1−bn<0 for all n≥n1. From (10), bn+1−bn=ηn+1(an+1−bn) and hence bn>an+1for n≥n1. Consequently, from (9), an+1>(1−ηn)2an+θ2an+1, which is equivalent to an+1>(1−ηn)2 1−θ2an. Since ηn→0, there is a µ > 1 and n2such that an+1≥µanfor all n≥n2. That is, (an)n≥n2grows at least geometrically, implying that after some time an+1>bn, which is a contradiction. To conclude, there is an m0≥n0such that (bn)n≥m0is strictly increasing. Lemma 10 shows that the expectation EuTSnuis ultimately bounded from below, assuming only that ηn→0. By additional assumptions on the sequence ηn, the growth rate can be characterised in terms of the adaptation weight sequence. 52
the random walk transition kernel with increment distribution q, and observe that the ‘adaptive random walk’ recursion (6) can be written as “Xn+1∼Q˜ qSn(Xn,·).” For any x∈Rdand measurable A⊂Rd |Q˜ qs(x,A)−P˜ qs(x,A)|≤21−Zmin1, π(y) π(x)˜ qs(y−x)dy ≤2 1−Z{|z|≤ ˜ M} min1, π(x+pθsz) π(x)˜ q(z)dz . Now, |Q˜ qs(x,A)−P˜ qs(x,A)| ≤ δwhenever pθsz ≤˜ δ1for all |z| ≤ ˜ M. In other words, there exists a µ=µ(δ)>0 such that whenever s< µ, the total variation norm kQ˜ qs(x,·)−P˜ qs(x,·)k≤δ. Next we shall consider a ‘adaptive random walk’ process to be coupled with (Xn,Mn,Sn)n≥1. Let n,k≥1 and define the random variables (˜ X(n) j,˜ M(n) j,˜ S(n) j)j∈[n,n+k]by setting (˜ X(n) n,˜ M(n) n,˜ S(n) n)≡ (Xn,Mn,Sn)and ˜ X(n) j+1∼Q˜ q˜ S(n) j (˜ X(n) j,·), ˜ M(n) j+1:= (1−ηj+1)˜ M(n) j+ηj+1˜ X(n) j+1and ˜ S(n) j+1:= (1−ηj+1)˜ S(n) j+ηj+1(˜ X(n) j+1−˜ M(n) j)2 for j+1∈[n+1,n+k]. The variable ˜ X(n) n+1can be selected so that P(˜ X(n) n+1=Xn+1| Fn) = 1−kP˜ qSn(Xn,·)−Q˜ q˜ S(n) n (˜ X(n) n,·)k; see Theorem 37 in Appendix D. Consequently, P(˜ X(n) n+16=Xn+1,Sn< µ|Fn)≤δ. By the same argument, ˜ X(n) n+2can be chosen so that P˜ X(n) n+26=Xn+2,˜ X(n) n+1=Xn+1,Sn+1< µ σ(Fn+1,˜ X(n) n+1)≤δ since if ˜ X(n) n+1=Xn+1, then also ˜ S(n) n+1=Sn+1. This implies P{˜ X(n) n+26=Xn+2}∪{˜ X(n) n+16=Xn+1}∩Bn:n+2Fn≤2δ where Bn:j:=∩j−1 i=n{Si< µ}for j>n. The same argument can be repeated to construct (˜ X(n) j)j∈[n,n+k]so that PDn:n+kFn≥1−kδ(17) where Dn:n+k:=Tn+k j=n{˜ X(n) j=Xj}∪B∁ n:n+k. Apply Lemma 16 with C=18 and ε=1/6 to obtain k≥1, and fix δ=ε/k. Denote ℓi:=ik +1 for any i≥0, and define the random variables (Ti)i≥1by Ti:= 1 {Sℓi−1<µ/2}min¨kMηℓi−1, ℓi X j=ℓi−1+1 logh1+ηjZ2 j−1i«(18) where Zjare defined as (13). 59
Define also ˜ Tisimilarly as Ti, but having ˜ Z(ℓi−1) jwith j∈[ℓi−1+1,ℓi]in the right hand side of (18), defined as ˜ Z(ℓi−1) ℓi−1≡Zℓi−1and by ˜ Z(ℓi−1) j:=˜ X(ℓi−1) j−˜ M(ℓi−1) j−1.q˜ S(ℓi−1) j−1. for j∈[ℓi−1+1,ℓi]. Notice that Ticoincides with ˜ Tiin Bℓi−1:ℓi∩Dℓi−1:ℓi. Observe also that ˜ X(ℓi−1) j follows the ‘adaptive random walk’ equation (6) for j∈[ℓi−1+1,ℓi], and hence ˜ Z(ℓi−1) jfollows (15). Consequently, denoting Gi:=Fℓi, Lemma 16 guarantees that PLℓi−1,kGi≤ε(19) where Lℓi−1,k:={˜ Ti<kMηℓi−1}. Let us show next that whenever Sℓi−1is small, the variable Tiis expected to have a positive value proportional to the adaptation weight, ETiGi−1 1 {Sℓi−1<µ/2}≥kηℓi−1 1 {Sℓi−1<µ/2}(20) almost surely for any sufficiently large i≥1. Write first ETiGi−1 1 {Sℓi−1<µ/2}=E( 1 B∁ ℓi−1:ℓi + 1 Bℓi−1:ℓi)TiGi−1 1 {Sℓi−1<µ/2} ≥E 1 B∁ ℓi−1:ℓi min§kCηℓi−1,µ 2+ξiª+ 1 Bℓi−1:ℓiξiGi−1 1 {Sℓi−1<µ/2} where the lower bound ξiof Tiis given as ξi:= ℓi X j=ℓi−1+1 log(1−ηj). By Assumption 15, ξi≥ −2kηℓi−1≥ −µ/4 for any sufficiently large i. Therefore, whenever PB∁ ℓi−1:ℓiGi−1≥ε=3/C, it holds that ETiGi−1 1 {Sℓi−1<µ/2}≥kηℓi−1 1 {Sℓi−1<µ/2} for any sufficiently large i. On the other hand, if PB∁ ℓi−1:ℓiGi−1≤ε, then by defining Ei:=B∁ ℓi−1:ℓi∪D∁ ℓi−1:ℓi∪Lℓi−1,k one has by (17) and (19) that P(Ei)≤3ε, and consequently ETiGi−1≥PE∁ iGi−1ξi+E 1 Ei˜ TiGi−1 ≥3εξi+ (1−3ε)kCηℓi−1≥kηℓi−1. 60
This establishes (20). Define the stopping times τ1≡1 and for n≥2 through τn:=inf{i> τn−1:Sℓi−1≥µ/2, Sℓi< µ/2} with the convention that inf;=∞. That is, τirecord the times when Sℓienters (0,µ/2]. Using τi, define the latest such time up to nby σn:=sup{τi:i≥1, τi≤n}. As in Theorem 18, define the almost surely converging martingale (Yi,Gi)i≥1with Y1≡0 and having the differences dYi:= (Ti−ETiGi−1)for i≥2. It is sufficient to show that liminfi→∞ Sℓi≥b:=µ/4>0 almost surely. If there is a finite i0≥1 such that Sℓi≥µ/2 for all i≥i0, the claim is trivial. Let us consider for the rest of the proof the case that {Sℓi< µ/2}happens for infinitely many indices i≥1. For any m≥2 such that Sℓm< µ/2, one can write logSℓm≥logSℓσm+ m X i=σm+1 Ti ≥logSℓσm+ (Ym−Yσm) + m X i=σm+1 kηℓi−1 (21) since then Sℓi< µ/2 for all i∈[σm,m−1]and hence also ETiGi−1≥kηℓi−1. Suppose for a moment that there is a positive probability that Sℓmstays within (0,µ/2)indefinitely, starting from some index m1≥1. Then, there is an infinite τiand consequently σm≤σ < ∞for all m≥1. But as Ymconverges, |Ym−Yσm|is a.s. finite, and since Pmηℓm=∞by Assumptions 15 and 17, the inequality (21) implies that Sℓm≥µ/2 for sufficiently large m, which is a contradiction. That is, the stopping times τifor all i≥1 must be a.s. finite, whenever Sℓm< µ/2 for infinitely many indices m≥1. For the rest of the proof, suppose Sℓm< µ/2 for infinitely many indices m≥1. Observe that since Ym→Y∞, there exists an a.s. finite index m2so that Ym−Y∞≥ −1/2log2 for all m≥m2. As ηn→0 and σm→ ∞, there is an a.s. finite m3such that ξσm−1≥ −1/2log2 for all m≥m3. For all m≥max{m2,m3}and whenever Sℓm< µ/2, it thereby holds that logSℓm≥logSℓσm−(Ym−Yσm)≥logSℓσm−1+ξσm−1 2log2 ≥log µ 2−log2 =log b. The case Sℓm≥µ/2 trivially satisfies the above estimate, concluding the proof. As a consequence of Theorem 19, one can establish a strong law of large numbers for the unconstrained AM algorithm running with a Laplace target distribution. Essentially, the only ingredient that needs to be checked is that the simultaneous geometric ergodicity condition holds. This is verified in the next lemma, whose proof is given in Appendix E. Lemma 20. Suppose that the template proposal distribution ˜ q is everywhere positive and nonincreasing away from the origin: ˜ q(z)≥˜ q(w)for all |z| ≤ |w|. Suppose also that π(x):= 1 2bexp−|x−m| bwith a mean m ∈Rand a scale b >0. Then, for all L >0, there are positive 61
constants M,b such that the following drift and minorisation condition are satisfied for all s ≥L and measurable A ⊂R PsV(x)≤λsV(x) + b 1 C(x),∀x∈R(22) Ps(x,A)≥δsν(A),∀x∈C(23) where V :R→[1,∞)is defined as V(x):= (supzπ(z))1/2π−1/2(x), the set C := [m−M,m+M], the probability measure µis concentrated on C and PsV(x):=RV(y)Ps(x,dy). Moreover, λs,δs∈(0,1) satisfy for all s ≥L max{(1−λs)−1,δ−1 s}≤ csγ(24) for some constants c,γ > 0that may depend on L. Theorem 21. Assume the adaptation weights (ηn)n≥2satisfy Assumptions 15 and 17, and the template proposal density ˜ q and the target distribution πsatisfy the assumptions in Lemma 20. If the functional f satisfies supx∈Rπ−γ(x)|f(x)|<∞for some γ∈(0,1/2). Then, n−1Pn k=1f(Xk)→Rf(x)π(x)dx almost surely as n →∞. Proof. The conditions of 19 are clearly satisfied implying that for any ε > 0 there is a κ=κ(ε)>0 such that the event Bκ:=§inf nSn≥κª has a probability P(Bκ)≥1−ε. The inequalities (22) and (23) of Lemma 20 with the bound (24) imply, using (Saksman and Vihola 2010, Proposition 7 and Lemma 12), that for any β > 0 there is a constant A=A(κ,ε,β)<∞such that P(Bκ∩{max{|Sn|,|Mn|}>Anβ})≤ε. Let us define the sequence of truncation sets Kn:={(m,s)∈R×R+:λmin(s)≥κ, max{|s|,|m|}≤Anβ} for n≥1. Construct an auxiliary truncated process (˜ Xn,˜ Mn,˜ Sn)n≥1, starting from (˜ X1,˜ M1,˜ S1)≡ (X1,M1,S1)and for n≥2 through ˜ Xn+1∼P˜ q˜ Sn(˜ Xn,·) (˜ Mn+1,˜ Sn+1) = σn+1h(˜ Mn,˜ Sn),ηn+1˜ Xn+1−˜ Mn,(˜ Xn+1−˜ Mn)2−˜ Sni where the truncation function σn+1:(Kn)×(R×R)→Knis defined as σn+1(z,z′) = (z+z′, if z+z′∈Kn z, otherwise. Observe that this constrained process coincides with the AM process with probability P∀n≥ 1 : (˜ Xn,˜ Mn,˜ Sn) = (Xn,Mn,Sn)≥1−2ε. Moreover, (Saksman and Vihola 2010, Theorem 2) implies that a strong law of large numbers holds for the truncated process (˜ Xn)n≥1, since supx|f(x)|V−α(x)<∞for some α∈(0,1 −β), by selecting β > 0 above sufficiently small. Since ε > 0 was arbitrary, the strong law of large numbers holds for (Xn)n≥1. 62
4 AM with a fixed proposal component This section deals with the modification due to Roberts and Rosenthal Roberts and Rosenthal (2009), including a fixed component in the proposal distribution. In terms of Section 2, the mixing parameter in (5) satisfies 0 < β < 1. Theorem 24 shows that the fixed proposal component guarantees, with a verifiable non-restrictive Assumption 22, that the eigenvalues of the adapted covariance parameter Snare bounded away from zero. As in Section 3.4, this result implies an ergodicity result, Theorem 29. Let us start by formulating the key assumption that, intuitively speaking, assures that the adaptive chain (Xn)n≥1will have ‘uniform mobility’ regardless of the adaptation parameter s∈Cd. Assumption 22. There exist a compactly supported probability measure νthat is absolutely continuous with respect to the Lebesgue measure, constants δ > 0 and c<∞and a measurable mapping ξ:Rd×Cd→Rdsuch that for all x∈Rdand s∈Cd, kξ(x,s)−xk≤ cand Pqs(x,A)≥δνA−ξ(x,s) for all measurable sets A⊂Rd, where A−y:={x−y:x∈A}is the translation of the set Aby y∈Rd. Remark 23.In the case of the AM algorithm with a fixed proposal component, one is primarily interested in the case where ξ(x,s) = ξ(x)and for all x∈Rd βqfix(x−y)min1, π(y) π(x)≥δνy−ξ(x) for all y∈Rd, where νis a uniform density on some ball. Then, since Pqs= (1−β)P˜ qs+βPqfix , Pqs(x,A)≥βPqfix (x,A)≥δZA ν(y−ξ)dy and Assumption 22 is fulfilled by the measure ν(A):=RAν(y)dy. Having Assumption 22, the lower bound on the eigenvalues of Sncan be obtained relatively easily, by a martingale argument similar to the one used in Section 3 and in Vihola (2009). Theorem 24. Let (Xn,Mn,Sn)n≥1be an AM process as defined in Section 2 satisfying Assumption 22. Moreover, suppose that the adaptation weights (ηn)n≥2satisfy Assumptions 15 and 17. Then, liminf n→∞ inf w∈S dwTSnw>0 where Sdstands for the unit sphere. Proof. Let us first introduce independent binary auxiliary variables (Zn)n≥2with Z1≡0, and through PZn+1=1Xn,Mn,Sn,Zn=δ PZn+1=0Xn,Mn,Sn,Zn= (1−δ). 63
Using this auxiliary variable, we can assume Xnto follow3 Xn+1=Zn+1(Un+1+ Ξn) + (1−Zn+1)Rn+1 where Un+1∼ν(·)is independent of Fnand Zn+1, the random variable Ξn:=ξ(Xn,Sn)is Fnmeasurable, and Rn+1is distributed according to the ‘residual’ transition kernel ˇ PSn(Xn,A):= (1− δ)−1[PqSn(Xn,A)−δν(A−Ξn)], valid by Assumption 22. Define S(w,γ):={v∈ S d:kw−vk ≤ γ}, the segment of the unit sphere centred at w∈ S dand having the radius γ > 0. Fix a unit vector w∈S dand define the following random variables Γ(γ) n+2:=inf v∈S (w,γ)|vT(Xn+1−Mn)|2+|vT(Xn+2−Mn+1)|2 for all n≥1. Denote Gn+1:=Xn+1−Mnand En+1:= Ξn+1−Xn+1, and observe that whenever Zn+2=1, it holds that Xn+2−Mn+1=Un+2+Xn+1−Mn+1+En+1 =Un+2+ (1−ηn+1)Gn+1+En+1 and we may write Zn+2Γ(γ) n+2=Zn+2inf v∈S (w,γ)|vTGn+1|2+|vT(Un+2+λn+1Gn+1+En+1)|2 where λn:=1−ηn∈(0,1)for all n≥2. Consequently, we may apply Lemma 25 below to find constants γ,µ > 0 such that PZn+2Γ(γ) n+2≥µFn≥δ 2. (25) Hereafter, assume γ > 0 is fixed such that (25) holds, and denote Γn+2:= Γ(γ) n+2and S(w):= S(w,γ). Consider the random variables Dn+2:=inf v∈S (w)ηn+1|vT(Xn+1−Mn)|2+ηn+2|vT(Xn+2−Mn+1)|2 ≥min{ηn+1,ηn+2}Γn+2≥η∗ηn+1Γn+2(26) where η∗:=infk≥2ηk+1/ηk>0 by Assumption 15. Define the indices ℓn:=2n−1 for n≥1 and let Tn:=η∗min{µ,ZℓnΓℓn} for all n≥2. Define the σ-algebras Gn:=Fℓnfor n≥1 and observe that ETn+1Gn≥ η∗µδ/2 by (25). Construct a martingale starting from Y1≡0 and having the differences dYn+1:= ηℓn+1(Tn+1−ETn+1Gn). The martingale Ynconverges to an a.s. finite limit Y∞as in Theorem 18. Define also η∗:=supk≥2ηk+1/ηk<∞and κ:=infk≥21−ηk>0, and let b:=κη∗µδ 8η∗>0. 3by possibly augmenting the probability space; see (Athreya and Ney 1978; Nummelin 1978). 64
Denote S(w) n:=infv∈S (w)vTSnvand define the stopping times τ1≡1 and for k≥2 through τk:=inf{n> τk−1:S(w) ℓn≤b,S(w) ℓn−1>b} with the convention inf;=∞. That is, τkrecord the times when S(w) ℓnenters (0, b]. Using τk, define the latest such time up to nby σn:=sup{τk:k≥1, τk≤n}. Observe that for any n≥2 such that S(w) ℓn≤b, one may write S(w) ℓn=S(w) ℓσn + n−1 X k=σnDℓk+2−ηℓk+1S(w) ℓk−ηℓk+2S(w) ℓk+1 ≥S(w) ℓσn + n−1 X k=σnηℓk+1Tk+1−ηℓk+1b−ηℓk+2κ−1b ≥S(w) ℓσn + n−1 X k=σn ηℓk+1Tk+1−η∗µδ 4 by (26) and since for all k∈[σn,n−1]one may estimate S(w) ℓk+1≤(1−ηℓk+1)−1S(w) ℓk+1≤κ−1b. That is, for any n≥2 such that S(w) ℓn≤b S(w) ℓn≥S(w) ℓσn + (Yn−Yσn) + n−1 X k=σn ηℓk+1ETk+1Gk−η∗δµ 4 ≥S(w) ℓσn + (Yn−Yσn) + η∗δµ 4 n−1 X k=σn ηℓk+1. As in the proof of Theorem 19, this is sufficient to find a ǫ > 0 such that liminf n→∞ S(w) n≥ǫ. Finally, take a finite number of unit vectors w1,...,wN∈S dsuch that the corresponding segments S(w1),...,S(wN)cover Sd. Then, liminf n→∞ inf v∈S dvTSnv=liminf n→∞ minS(w1) n,...,S(wN) n≥ǫ. Lemma 25. Suppose Fn⊂ Fn+1are σ-algebras, and Gn+1and En+1are Fn+1-measurable random variables, satisfying kEn+1k ≤ M for some constant M <∞. Moreover, Un+2is a random variable independent of Fn+1, having a distribution νfulfilling the conditions in Assumption 22. Let Sd:={u∈Rd:kuk=1}stand for the unit sphere and denote by S(w,γ):={v∈S d:kw−vk≤ γ}the segment of the unit sphere centred at w ∈S dand having the radius γ > 0. There exist constants γ,µ > 0such that Pinf v∈S (w,γ)|vTGn+1|2+|vT(Un+2+λGn+1+En+1)|2> µ Fn≥1 2. for any w ∈S dand any constant λ∈(0,1), almost surely. 65
Proof. Since νis absolutely continuous with respect to the Lebesgue measure, one can show that there exist values b,γ > 0 such that inf w∈S dinf e∈B(0,M)νu∈Rd: inf v∈S (w,γ)|vT(u+e)|>b≥1 2(27) where B(0, M):={y∈Rd:kyk ≤ M}denotes a centred ball of radius M. Hereafter, fix γ,b>0 such that (27) holds and let a:=b/2. Fix a unit vector w∈S dand consider the set A:=inf v∈S (w,γ)|vTGn+1|2+|vT(Un+2+λGn+1+En+1)|2≤a2 ⊂¨inf v∈S (w,γ):|vTGn+1|≤a|vT(Un+2+λGn+1+En+1)|≤ a« ⊂¨inf v∈S (w,γ):|vTGn+1|≤a|vT(Un+2+En+1)|−λ|vTGn+1|≤ a« ⊂inf v∈S (w,γ)|vT(Un+2+En+1)|≤ 2a. Since Un+2is independent of Fn+1, and since En+1is Fn+1-measurable, one may estimate PA∁Fn≥E inf e∈B(0,M) Pinf v∈S (w,γ)|vT(Un+2+e)|>2aFn+1Fn =inf e∈B(0,M)νu∈Rd: inf v∈S (w,γ)|vT(u+e)|>b≥1 2 by (27), almost surely, concluding the proof by µ:=a2. Corollary 26. Assume πis bounded, stays bounded away from zero on compact sets, is differentiable on the tails, and has regular contours, that is, liminf kxk→∞ x kxk·∇π(x) k∇π(x)k<0. (28) Let (Xn,Mn,Sn)n≥1be an AM process as defined in Section 2 using a mixture proposal (5) with a mixing weight satisfying β∈(0,1)and the density qfix is bounded away from zero in some neighbourhood of the origin. Moreover, suppose that the adaptation weights (ηn)n≥2satisfy Assumptions 15 and 17. Then, liminf n→∞ inf w∈S dwTSnw>0. Proof. In light of Theorem 24, it is sufficient to check Assumption 22, or in fact the conditions in Remark 23. Let L>0 be sufficiently large so that infkxk≥Lx kxk·∇π(x) k∇π(x)k<0. Jarner and Hansen (Jarner and Hansen 2000, proof of Theorem 4.3) show that there is an ε′>0 and K>0 such that the cone E(x):=¨x−au : 0 <a<K,u∈S d, u−x kxk ≤ε′« 66
is contained in the set A(x):={y∈Rd:π(y)≥π(x)}, for all kxk≥ L. Let r′>0 be sufficiently small to ensure that infkzk≤r′qfix(z)≥δ′>0. There is a r=r(ε′,K)∈ (0, r′/2)and measurable ξ:Rd→Rdsuch that kξ(x)−xk ≤ r′/2 and the ball B(x,r):={y: ky−ξ(x)k ≤ r}is contained in the cone E(x). Define ν(x):=c−1 r 1 B(0,r)(x)where cr:=|B(0, r)| is the Lebesgue measure of B(0, r), and let ξ(x):=xfor the remaining kxk<L. Now, we have for kxk≥ Lthat βqfix(x−y)min1, π(y) π(x)≥βδ′crν(y−ξ). Since πis bounded and bounded away from zero on compact sets, the ratio π(y)/π(x)≥δ′′ >0 for all x,y∈B(0, L+r′)with kx−yk≤ r′. Therefore, for all kxk<L, it holds that βqfix(x−y)min1, π(y) π(x)≥βδ′δ′′crν(y−x). Remark 27.The conditions of Corollary 26 are fulfilled by many practical densities π(see Jarner and Hansen (2000) for examples), and are fairly easy to verify in practice. Assumption 22 holds, however, more generally, excluding only densities with unbounded density or having irregular contours. Remark 28.It is not necessary for Theorem 24 and Corollary 26 to hold that the adaptive proposal densities {˜ qs}s∈Cdhave the specific form discussed in Section 2. The results require only that a suitable fixed proposal component is used so that Assumption 22 holds. In Theorem 29 below, however, the structure of {˜ qs}s∈Cdis required. Let us record the following ergodicity result, which is a counterpart to (Saksman and Vihola 2010, Theorem 10) formulating a a strong law of large numbers for the original algorithm (S1)–(S3) with the covariance parameter (1). Theorem 29. Suppose the target density πis continuous and differentiable, stays bounded away from zero on compact sets and has super-exponentially decaying tails with regular contours, limsup kxk→∞ x kxkρ·∇logπ(x) = −∞ and limsup kxk→∞ x kxk·∇π(x) k∇π(x)k<0, respectively, for some ρ > 1. Let (Xn,Mn,Sn)n≥1be an AM process as defined in Section 2 using a mixture proposal qs(z) = (1− β)˜ qs(z) + βqfix(z)where ˜ qsstands for a zero-mean Gaussian density with covariance s, the mixing weight satisfies β∈(0,1)and the density qfix is bounded away from zero in some neighbourhood of the origin. Moreover, suppose that the adaptation weights (ηn)n≥2satisfy Assumption 17. Then, for any function f :Rd→Rwith supx∈Rdπγ(x)|f(x)|<∞for some γ∈(0,1/2), 1 n n X k=1 f(Xk)n→∞ −−−→ZRd f(x)π(x)dx almost surely. 67
Proof. The conditions of Corollary 26 are satisfied, implying that for any ε > 0 there is a κ=κ(ε)> 0 such that Pinfnλmin(Sn)≥κ≥1−εwhere λmin(s)denotes the smallest eigenvalue of s. By (Saksman and Vihola 2010, Proposition 15), there is a compact set Cκ⊂Rd, a probability measure νκon Cκ, and bκ<∞such that for all s∈ Cdwith λmin(s)≥κ, it holds that P˜ qsV(x)≤λsV(x) + b 1 Cκ(x),∀x∈Rd(29) P˜ qs(x,A)≥δsν(A)∀x∈Cκ(30) where V(x):= (supxπ(x))1/2π−1/2(x)≥1 and the constants λs,δs∈(0,1)satisfy the bound (1−λs)−1∨δ−1 s≤c1det(s)1/2(31) for some constant c1≥1. Likewise, there is a compact Df⊂Rd, a probability measure µfon Df, and constants bf<∞and λf,δf∈(0,1), so that (29) and (30) hold with Pf(Jarner and Hansen 2000, Theorem 4.3). Put together, (29) and (30) hold for Pqsfor all s∈ Cdwith λmin(s)≥κ, perhaps with different constants, but satisfying a bound (31), with another c2≥c1. The rest of the proof follows as in Theorem 21 by construction of an auxiliary process (˜ Xn,˜ Mn,˜ Sn)n≥1 truncated so that for given ǫ > 0, κ≤λmin(˜ Sn)≤anǫand |˜ Mn| ≤ anǫand where the constant a=a(ǫ,κ)is chosen so that the truncated process coincides with the original AM process with probability ≥1−2ε. Theorem 2 of Saksman and Vihola (2010) ensures that the strong law of large numbers holds for the constrained process, and letting ε→0 implies the claim. Remark 30.In the case ηn:=n−1, Theorem 29 implies that with probability one, Mn→mπ:= Rxπ(x)dxand Sn→sπ:=Rx xTπ(x)dx−mπmT π, the true mean and covariance of π, respectively. Remark 31.Theorem 29 holds also when using multivariate Student distributions {˜ qs}s∈Cd, as (Vihola 2009, Proposition 26 and Lemma 28) extend the result in Saksman and Vihola (2010) to cover this case. Acknowledgements The author thanks Professor Eero Saksman for discussions and helpful comments on the manuscript. References C. Andrieu and É. Moulines. On the ergodicity properties of some adaptive MCMC algorithms. Ann. Appl. Probab., 16(3):1462–1505, 2006. MR2260070 C. Andrieu and C. P. Robert. Controlled MCMC for optimal sampling. Technical Report Ceremade 0125, Université Paris Dauphine, 2001. C. Andrieu and J. Thoms. A tutorial on adaptive MCMC. Statist. Comput., 18(4):343–373, Dec. 2008. MR2461882 Y. Atchadé and G. Fort. Limit theorems for some adaptive MCMC algorithms with subgeometric kernels. Bernoulli, 16(1):116–154, Feb. 2010. MR2648752 68
where a(x,y):= 1−Èπ(x) π(y) =1−e−x−|y| 2and b(x,y):=rπ(y) π(x) 1−rπ(y) π(x)!=e−|y|−x 21−e−|y|−x 2. Compute then that Zx 0 a(x,y)˜ qs(y−x)dy−Z2x x b(x,y)˜ qs(y−x)dy=Zx 01−e−z 22˜ qs(z)dz. The estimates Z0 −x a(x,y)˜ qs(y−x)dy≥˜ qs(2x)Zx 0 a(x,y)dy=˜ qs(2x)Zx 0 (1−e−z 2)dz Z−x −∞ b(x,y)˜ qs(y−x)dy≤˜ qs(2x)Z∞ x b(x,y)dy=˜ qs(2x)Z∞ 0 e−z 2(1−ez 2)dz due to the non-increasing property of ˜ qsyield Z0 −x a(x,y)˜ qs(y−x)dy−Z−x −∞ b(x,y)˜ qs(y−x)dy ≥˜ qs(2x)Zx 0 (1−e−z 2)2dz−Z∞ x e−z 2dz>0 for any sufficiently large x>0. Similarly, one obtains 1 2Zx 01−e−z 22˜ qs(z)dz−Z∞ 2x b(x,y)qs(y−x)dy>0 for large enough x>0. Summing up, letting M>0 be sufficiently large, then for x≥Mand s≥L>0 1−PsV(x) V(x)≥1 2Zx 01−e−z 22˜ qs(z)dz≥1 2˜ qs(M)ZM 01−e−z 22dz ≥c1s−1/2˜ q(θ−1/2s−1/2M)≥c2s−1/2 for some constants c1,c2>0. The same inequality holds also for −x≤ −Mdue to symmetry. The simple bound PsV(x)≤2V(x)observed from (32) with the above estimate establishes (22). The minorisation inequality (23) holds since for all x∈Cone may write Ps(x,A)≥ZA∩C max1, π(y) π(x)˜ qs(y−x)dy ≥infz∈Cπ(z) supzπ(z)inf s≥L,z,y∈C˜ qs(z−y)ZA∩C dy≥c3s−1/2ν(A). where ν(A):=|A∩C|/|C|with |·| denoting the Lebesgue measure. 75