Efficiency and computability of MCMC with Langevin, Hamiltonian, and other matrix-splitting proposals

0.31 MB·English

by  Richard A. Norton

Arxivversionpreprint
EFFICIENCY AND COMPUTABILITY OF MCMC WITH
LANGEVIN, HAMILTONIAN, AND OTHER
MATRIX-SPLITTING PROPOSALS
5
By Richard A. Norton and Colin Fox
1
0
University of Otago
2

n WeanalysecomputationalefficiencyofMetropolis-Hastingsalgo-
a rithms with AR(1) process proposals. These proposals include, as a
J subclass,discretized Langevindiffusion(e.g.MALA)anddiscretized

3 Hamiltonian dynamics(e.g. HMC).
1 By including the effect of Metropolis-Hastings we extend earlier
work by Fox and Parker, who used matrix splitting techniques to
]
analyse the performance and improve efficiency of AR(1) processes
R
for targeting Gaussian distributions.
P
Ourresearch enablesanalysisof MCMC methodsthatdrawsam-
.
h ples from non-Gaussian target distributions by using AR(1) process
t proposalsinMetropolis-Hastingsalgorithms,byanalysingthematrix
a
m splitting of the precision matrix for a local Gaussian approximation
of the non-Gaussian target.
[

1
v 1. Introduction. Many existingMetropolis-Hastings (MH)algorithms
0 forsamplingfromanon-Gaussiantargetdistributionπ useanAR(1)process
5
proposal; given current state x Rd the proposal y Rd is given by
1
∈ ∈
3
0 (1) y = Gx+g+ν
.
1
0 where G Rd×d is the iteration matrix, g Rd is a fixed vector and ν is an
∈ ∈
5 i.i.d. draw from N(0,Σ). In general, G, g and Σ may depend on x. The MH
1
acceptance probability is
:
v
i π(y)q(y,x)
X α(x,y) = 1
∧ π(x)q(x,y)
r
a
where π(x) denotes the target probability density function and q(x,dy) =
q(x,y)dy is the transition kernel for the proposal y given current state x.
An important case of (1) is when the AR(1) process comes from a matrix
splittingofsomematrix associated withπ atx,forexamplewhen isthe
A A
MSC 2010 subject classifications: Primary 60J22, 65Y20, ; secondary 60J05, 62M05,
65C40, 68Q87, 68W20
Keywords and phrases: Markov chain, Monte Carlo, sampling, Gaussian, multivariate
normal distribution, AR(1) process, Metropolis-Hastings algorithm, Metropolis-adjusted
Langevin algorithm, Hybrid Monte Carlo,
1
2 R. A.NORTONAND C. FOX
local Hessian of logπ, or its inverse, or an approximation to these matrices.
In this case we find it more convenient and revealing to write (1) in terms
of the matrix splitting = M N so that y is given by solving
A −
(2) My = Nx+β +ν
where β is a vector and ν is an i.i.d. draw from N(0,MT +N), G =M 1N,
−
g = M 1β and Σ = M 1(MT +N)M T. In fact, if the spectral radius of
− − −
G is less than 1 then the AR(1) processes (1) can be written in terms of
a matrix splitting and vice versa, see Section 2. Moreover, if the spectral
radius of G is less than 1, then the proposal chain generated by (2) will
converge to N( 1β, 1), which we call the proposal target distribution,
− −
A A
see [12].
Later, we see that the Metropolis-adjusted Langevin algorithm (MALA)
is an example of an AR(1) process proposal with a corresponding matrix
splitting,asareotherdiscretizedLangevindiffusionproposals,asarealsothe
Hybrid Monte Carlo algorithm (HMC) and other discretized Hamiltonian
dynamics proposals. We will use this observation to further analyse these
methods and discuss their computability in Section 5.
Although identifying (1)with (2)is usefulforanalysing existingmethods,
if the task is designing a proposal for the MH algorithm then using (2) is
morenaturalbecausewecanbeginbychoosing andβ.Forexample,weare
A
particularly interested in the proposal where and β are chosen such that
A
N( 1β, 1) is a local Gaussian approximation to π, i.e. 1xT x+βTx
A− A− −2 A
is a local quadratic approximation to logπ.
After selecting and β,we mustalsochoose the splittingM and N.This
A
choice willeffect thecomputational costof each proposal,as wellashow fast
the AR(1) process will converge to the proposal target distribution.
We may use one or more iterations of the AR(1) process (2) as a pro-
posal for the MH algorithm. For the next proposal (if the previous one was
accepted) we then use a new local Gaussian approximation to π. This idea
mimics thedesignofsomeoptimizers wherealocal quadraticapproximation
to theobjective functionisminimized ateach iteration, seeforexample[18].
Ourconcern is the analysis of Markov chain Monte Carlo(MCMC) meth-
ods with AR(1) process proposals which fall into three cases.
1. The target and proposal target distributions are the same, i.e. π =
N( 1β, 1) and the AR(1) process targets π directly. Then (1) is
− −
A A
a fixed-scan Gibbs sampler for π, and the accept/reject step in the
MH algorithm is redundant. Fox and Parker [11, 12, 13] have already
studied this case. They showed how to choose M and N (equivalently
EFFICIENCY OF MCMC 3
choose G and Σ) to obtain an efficient algorithm and they also showed
how to accelerate the AR(1) process to achieve even greater efficiency.
2. The target distribution π is normal, N(A 1b,A 1) for some matrix
− −
A and vector b, but not targeted by the AR(1) process, i.e. A =
6 A
and/or b = β. Then a MH accept/reject step is required to ensure π
6
is targeted.
3. The target distribution π is not normal, hence the AR(1) process does
not target π, and a MH accept/reject step is required to ensure π is
targeted. The proposal target is a local Gaussian approximation to π.
Cases 2 and 3 where the target and proposal target distributions are
different are harder to analyse than case 1 as the MH accept/reject step
affects thetransition kernel.Thisbegsthefollowingquestions inthesecases:
How should we choose , β, M and N to construct a computationally
• A
efficient algorithm?
Can we accelerate MH algorithms with AR(1) process proposals for
•
non-Gaussian target distributions?
We measurethequality of asamplingalgorithm byits computational effi-
ciency, which is a combination of the compute time required for the Markov
chain to reach equilibrium (burn in), and once in equilibrium, the compute
time required to produce quasi-independent samples. Simply, we are inter-
ested in getting the most ‘bang for buck’, where ‘bang’ is the number of
quasi-independent samples and ‘buck’ is the computing cost.
Fox and Parker [12] mainly considered distributional convergence from
an arbitrary state, so they measured burn in. Here we are concerned about
the compute time to produce quasi-independent samples. The integrated
autocorrelation time (length of Markov chain with variance reducing power
of one independent sample) is a proxy for independence, hence we try to
measure compute time per integrated autocorrelation time.
In this article we focus on analysing case 2 above where π is the nor-
mal N(A 1b,A 1), but not the same as the proposal target N( 1β, 1).
− − − −
A A
However, our results are relevant to case 3 where π is not normal because
the cases share important features. In particular, the accept/reject step in
the MHalgorithm is usedin both cases to correct thedifference between the
target andproposaltarget distributions.Thesimplification ofonly consider-
ingcase2makes theanalysisofthetransition kernelpossible,whilstkeeping
the essential feature of the non-Gaussian target case, that the spectrum of
the transition kernel has been changed by the MH algorithm. If we cannot
analyse or accelerate case 2 then we do not expect to be able to analyse or
accelerate case 3.
4 R. A.NORTONAND C. FOX
Moreover, in case 3 we might construct an inhomogeneous Markov chain
by updating the proposal rule at each iteration (or every few iterations) by
updating the local Gaussian approximation of the target distribution, hence
we are interested in the local behaviour of the algorithm where the target
distribution is approximately normal.
Our analysis of case 2 is for the special case when the matrices M and
N are functions of A. This allows us to simultaneously diagonalise both the
AR(1) process and the target distribution with a change of variables. This
is not an overly restrictive condition if factorizing the matrix A is infeasible,
and we will see below that it includes several important examples of MH
algorithms already in use.
Our analysis is also analogous to the analysis of optimizers where it is
useful to test the behaviour of optimizers on quadratic cost functions (in
the probability sense, a Gaussian is equivalent to a quadratic function since
logπ is quadratic when π is Gaussian).
Acceleration techniques for AR(1) processes for sampling in case 1 (nor-
mal target) do not necessarily accelerate sampling in case 3 (non-normal
target) when the accelerated AR(1) process is used as a proposal in the
MH algorithm. Goodman and Sokal [14] accelerated Gibbs sampling of nor-
mal distributions using ideas from multigrid linear solvers, but only ob-
served modest efficiency gains in their non-normal examples (exponential
distributions with fourth moments). Green and Han [15] applied successive-
over-relaxation to a local Gaussian approximation of a non-Gaussian target
distributionasaproposalfortheMHalgorithm.Again,theydidnotobserve
significant acceleration in the non-normal target case.
By recognising that MH algorithms with discretized Langevin diffusion
or discretized Hamiltonian dynamics proposals are examples of MH algo-
rithms with AR(1) process proposals, we are able to apply our theory to
these MCMC methods. We find that we can replicate existing results in the
literature for these methods,and we can also extend existing results in some
cases, see Section 5. We also discuss the computational efficiency of these
methods.
Strictly, our analysis is exact for Gaussian target distributions, but it can
also be applied to the case when π is absolutely continuous with respect to
a Gaussian (see e.g. [24], [5] or [6]).
Sampling from a Gaussian (or a target that is absolutely continuous with
respecttoaGaussian)isalsoaproblemthatarisesininverseproblemswhere
the forward model has the form y = Fx+n where y is the observed data, F
isalinear(orboundedforabsolutely continuous withrespecttoaGaussian)
operator, x is the unknown image with Gaussian prior distribution, and n is
EFFICIENCY OF MCMC 5
Gaussian noise, see e.g. [24, Example 6.23]. If the posterior is Gaussian with
hyperparameters satisfying some other marginal distribution, and we can
sample the hyperparameters independently,then conditionally sampling the
posterior given the hyperparameters is also a Gaussian sampling problem
[2]. Other applications that involve sampling Gaussian distributions include
Brownian bridges [6].
For the purpose of evaluating the cost of computing a proposal from (2)
for agiven matrix splitting,wemake similarassumptionsas madein numer-
ical linear algebra for solving Ax = b, since solver and sampling algorithms
share many of the same operations. We assume that we can efficiently com-
pute matrix-vector products Av for any vector v Rd; for example A may
∈
be sparse. Furthermore, we assume that d is sufficiently large so that it is
infeasible to directly compute A1/2 or any other matrix factorization of A
and then directly compute independent samples from our Gaussian target.
The remaining sections are as follows. In Section 2 we quickly show that
(1) and (2) are equivalent, then Sections 3 and 4 present new analyses for
the expected acceptance rate and jump size of MH algorithms with AR(1)
processproposals.Section 5then appliesthis newanalysisto proposalsfrom
Langevin diffusion and Hamiltonian dynamics. We see that these proposals
are AR(1) processes and we identify the corresponding matrix splitting and
proposal target distribution. Using our earlier analysis we assess the con-
vergence properties of these methods as d . We provide concluding
→ ∞
remarks in Section 6.
2. AR(1) processes correspond to matrix splittings. We can ex-
press an AR(1) process using either (1) or (2), provided the AR(1) process
converges.
Theorem 2.1. If we are given G, g and Σ, and the spectral radius of G
is less than 1, then the AR(1) process (1) can be written as (2) using
1
∞ −
(3) = GlΣ(GT)l
A
!
l=0
X
and
M = (I G) 1,
−
(4) A −
N = (I G) 1G, β = (I G) 1g.
− −
A − A −
Note that we also have = M N symmetric and positive definite.
A −
6 R. A.NORTONAND C. FOX
Proof. Since the spectral radius of G is less than 1 and Σ is symmet-
ric positive definite, it follows that A−1 := ∞l=0GlΣ(GT)l is well-defined,
symmetric and positive definite.
P
Itis theneasy toseethat(3)and(4)satisfy = M N,G = M 1N and
−
A −
g = M 1β. We must also check that Σ = M 1(MT +N)M T is satisfied.
− − −
Substituting (4) into M 1(MT +N)M T we get
− −
M 1(MT +N)M T = (I G) 1+G 1GT
− − − −
− A A
= 1 G 1GT
− −
A − A
∞
= GlΣ(GT)l GlΣ(GT)l
−
l=0 l=1
X X
= Σ.
In the following special case we obtain a symmetric matrix splitting.
Corollary 2.2. If the spectral radius of G is less than 1 and GΣ is
symmetric, then the AR(1) process (1) has a corresponding matrix splitting
defined by
M = Σ 1(I +G), = M(I G) = Σ 1(I G2),
− −
A − −
N = MG = Σ 1(I +G)G, β = Mg = Σ 1(I +G)g,
− −
and M and N are symmetric (we say the matrix splitting is symmetric).
Proof. These matrix splitting formulae follow from the identity
∞ ∞
GlΣ(GT)l = G2lΣ = (I G2) 1Σ.
−
−
l=0 l=0
X X
To see that M is symmetric (and hence also N since is symmetric and
A
= M N) we note that
A −
1
∞ −
M = (I G) G2lΣ .
−
!
l=0
X
EFFICIENCY OF MCMC 7
3. Expected acceptance rate for a general matrix splitting. The
expected acceptance rate is a quantity that is related to efficiency and for
optimal performance the proposal is usually tuned so that the observed
average acceptance rate is well away from 0 and 1. For example, it has been
shown in the case when d that 0.234 is optimal for the random-walk
→ ∞
Metropolis algorithm (RWM) [20], 0.574 is optimal for MALA [21], and
0.651 is optimal for HMC [3]. All of these results required expressions for
the expected acceptance rate of the algorithm as d . Here we derive an
→ ∞
expression for the expected acceptance rate for an AR(1) process proposal
(2) with Gaussian target N(A 1b,A 1), provided the splitting matrices are
− −
functions of A. Thus, our MH algorithm is defined by
Target: N(A 1b,A 1),
− −
(5)
Proposal: y = Gx+M 1β +( 1 G 1GT)1/2ξ,
− − −
A − A
whereξ N(0,I),andwehaveusedG = M 1N andM 1(MT+N)M T =
− − −
∼
1 G 1GT [13, Lem. 2.3]. The following lemma is a result of simple
− −
A − A
algebra. The proof is in the Appendix.
Lemma 3.1. Suppose = M N is a symmetric splitting. Then the
A −
MH acceptance probability for (5) satisfies
1 1
α(x,y) = 1 exp yT(A )y+ xT(A )x+(b β)T(y x) .
∧ −2 −A 2 −A − −
(cid:18) (cid:19)
Now define a spectral decomposition
A= QΛQT
whereQ Rd d is orthogonal and Λ= diag(λ2,...,λ2) is a diagonal matrix
∈ × 1 d
of eigenvalues of A. Although we may not be able to compute Q and Λ this
does not stop us from using the spectral decomposition for theory. Simple
algebra gives us the following result.
Lemma 3.2. Suppose M = M(A) and N = N(A) are functions of A.
Then G and are also functions of A and under the coordinate transfor-
A
mation
x QTx
↔
the MH algorithm (5) is transformed to a MH algorithm defined by
Target: N(Λ 1QTb,Λ 1),
− −
(6)
Proposal: y = Gx+M(Λ) 1QTβ+(Λ˜ 1 GΛ˜ 1GT)1/2ξ,
− − −
−
8 R. A.NORTONAND C. FOX
where ξ N(0,I), and G = M(Λ) 1N(Λ) and Λ˜ = (Λ) are diagonal
−
∼ A
matrices.
UsingLemma3.1weseethattheacceptanceprobabilityofMHalgorithms
(5) and (6) are identical and hence it is sufficient to analyse the convergence
properties of (6) to determine the convergence properties of (5).
We will need the following Lyapunov version of the Central Limit Theo-
rem, see e.g. [7, Thm. 27.3]
Theorem 3.3. Let X ,...,X be a sequence of independent random
1 d
variables withfiniteexpected valueµ and variance σ2. Defines2 := d σ2.
i i d i=1 i
If there exists a δ > 0 such that
P
d
1
lim E[X µ 2+δ] = 0,
d s2+δ | i − i|
→∞ d i=1
X
then
d
1
(Xi µi) D N(0,1) as d .
s − −→ → ∞
d
i=1
X
An equivalent conclusion to this theorem is d X N( d µ ,s2) in
i=1 i → i=1 i d
distribution as d . Another useful fact is
→ ∞ P P
(7) X N(µ,σ2) E[1 eX] = Φ(µ)+eµ+σ2/2Φ( σ µ)
∼ ⇒ ∧ σ − − σ
where Φ is the standard normal cumulative distribution function. See e.g.
[20, Prop. 2.4] or [5, Lem. B.2].
Theorem 3.4. Suppose that M and N in (5) are functions of A, and
the Markov chain is in equilibrium, i.e. x N(A 1b,A 1). Recall that λ2
are eigenvalues of A, and define λ˜2 and G∼to be t−he eig−envalues of andi
i i A
G respectively. Also define
g˜ := 1 G , g := 1 G2, m := (A 1b) , m˜ := ( 1β) ,
i − i i − i i − i i A− i
λ2 λ˜2 λ2 λ˜2
r := i − i , r˜ := i − i , rˆ := m m˜
i λ2 i λ˜2 i i − i
i i
EFFICIENCY OF MCMC 9
and
T := rˆ2λ2(1r g g˜),
0i i i 2 i i − i
T := rˆλ (r g g˜),
1i i i i i i
−
T := rˆλ g1/2(1+r˜)1/2(1 r G ),
2i i i i i − i i
T := 1r g ,
3i 2 i i
T := 1r g (1+r˜),
4i −2 i i i
T := r G g1/2(1+r˜)1/2.
5i − i i i i
If there exists a δ > 0 such that
d T 2+δ
(8) lim i=1| ji| = 0 for j = 1,2,3,4,5
d→∞ Pdi=1|Tji|2 1+δ/2
(cid:16) (cid:17)
(j = 0 is not requirPed), then
π(y)q(y,x)
Z := log D N(µ,σ2) as d
π(x)q(x,y) −→ → ∞
(cid:18) (cid:19)
where µ = lim d µ and σ2 = lim d σ2 and
d i=1 i d i=1 i
→∞ →∞
µi = T0i +T3i +PT4i and σi2 = T12i +PT22i+2T32i +2T42i +T52i.
Hence, the expected acceptance rate satisfies
E[α(x,y)] Φ(µ)+eµ+σ2/2Φ( σ µ) as d
→ σ − − σ → ∞
where Φ is the standard normal cumulative distribution function.
Proof. By Lemma 3.2 it is sufficient to only consider (5) in the case
where all matrices are diagonal matrices, e.g. A = diag(λ2,...,λ2), =
diag(λ˜2,...,λ˜2),M = diag(M ,...,M ),G = diag(G ,...,G1 ),m d=λA2b ,
1 d 1 d 1 d i −i i
and m˜ = λ˜ 2β . Then, in equilibrium we have
i −i i
1
x = m + ξ
i i i
λ
i
where ξ N(0,1) and using m˜ =G m˜ +M 1β we have
i ∼ i i i i− i
β (1 G2)1/2
y = G x + i + − i ν
i i i Mi λ˜i i
1/2
1 g
= G m + ξ +(1 G )m˜ + i ν
i(cid:18) i λi i(cid:19) − i i λ˜i i
1/2
G g
= m˜ +G rˆ+ iξ + i ν
i i λi i λ˜i i
10 R. A.NORTONAND C. FOX
where ν N(0,1). From Lemma 3.1 we also have Z = d Z where
i ∼ i=1 i
Z = 1(λ2 λ˜2)(y2 x2)+(b β )(y Px ).
i −2 i − i i − i i − i i− i
Substituting x and y as above, using identities such as λ2λ˜ 2 =1+r˜ and
i i i −i i
(b β )λ 2 =rˆ +r m˜ , then after some algebra we eventually find
i − i −i i i i
Z =T +T ξ +T ν +T ξ2+T ν2+T ξ ν .
i 0i 1i i 2i i 3i i 4i i 5i i i
Hence
µ := E[Z ]= T +T +T
i i 0i 3i 4i
and
σ2 := Var[Z ]= E[Z2] E[Z ]2
i i i − i
= T2 +T2 +T2 +3T2 +3T2 +T2 +2T T +2T T +2T T
0i 1i 2i 3i 4i 5i 0i 3i 0i 4i 3i 4i
(cid:0) (T0i +T3i +T4i)2 (cid:1)
−
= T2 +T2 +2T2 +2T2 +T2
1i 2i 3i 4i 5i
and
Z µ = T ξ +T ν +T (ξ2 1)+T (ν2 1)+T ξ ν .
i− i 1i i 2i i 3i i − 4i i − 5i i i
Therefore, for any d N and δ > 0 we can bound the Lyapunov condition
∈
in Theorem 3.3 as follows
1 d 52+δ 5 d
E[Q µ 2+δ] C (δ) T 2+δ
s2+δ | i − i| ≤ s2+δ j | ji|
d i=1 d j=1 i=1
X X X
5 d T 2+δ
52+δ C (δ) i=1| ji|
j
≤ 1+δ/2
j=1 Pd T2
X i=1 ji
(cid:16) (cid:17)
P
where C (δ) = C (δ) = E[ξ 2+δ] , C (δ) = C (δ) = E[ξ2 12+δ] and
1 2 3 4
| | | − |
C (δ) =E[ξ 2+δ]2, and ξ N(0,1).
5
| | ∼
Therefore, if (8) holds then the result follows from Theorem 3.3 and (7).
4. Expected squared jump size for a general matrix splitting.
Theefficiency ofaMCMCmethodisusuallygivenintermsoftheintegrated
autocorrelation time which may be thought of as “the number of dependent
sample points from the Markov chain needed to give the variance reducing
power of one independent point”, see e.g. [17, 6.3]. Unfortunately, we are
§