Estimating the spectral gap of a trace-class Markov operator

arXiv:1704.00850v2 [] 22 Jun 2017

Qian Qin, James P. Hobert and Kshitij Khare
Department of Statistics
University of Florida
June 2017

Abstract
The utility of a Markov chain Monte Carlo algorithm is, in large part, determined by the size of
the spectral gap of the corresponding Markov operator. However, calculating (and even approximating)
the spectral gaps of practical Monte Carlo Markov chains in statistics has proven to be an extremely
difficult and often insurmountable task, especially when these chains move on continuous state spaces.
In this paper, a method for accurate estimation of the spectral gap is developed for general state space
Markov chains whose operators are non-negative and trace-class. The method is based on the fact that the
second largest eigenvalue (and hence the spectral gap) of such operators can be bounded above and below
by simple functions of the power sums of the eigenvalues. These power sums often have nice integral
representations. A classical Monte Carlo method is proposed to estimate these integrals, and a simple
sufficient condition for finite variance is provided. This leads to asymptotically valid confidence intervals
for the second largest eigenvalue (and the spectral gap) of the Markov operator. The efficiency of the
method is studied. For illustration, the method is applied to Albert and Chib’s (1993) data augmentation
(DA) algorithm for Bayesian probit regression, and also to a DA algorithm for Bayesian linear regression
with non-Gaussian errors (Liu, 1996).

1 Introduction
Markov chain Monte Carlo (MCMC) is widely used to estimate intractable integrals that represent expectations with respect to complicated probability distributions. Let π : S → [0, ∞) be a probability density
function (pdf) with respect to a σ-finite measure µ, where (S, U , µ) is some measure space. Suppose we
Key words and phrases. Data augmentation algorithm, Eigenvalues, Hilbert-Schmidt operator, Markov chain, Monte Carlo

1

want to approximate the integral
J :=

Z

f (u)π(u)µ(du)

S

P
m−1
for some function f : S → R. Then J can be estimated by Jˆm := m−1 m−1
k=0 f (Φk ), where {Φk }k=0 are

the first m elements of a well-behaved Markov chain with stationary density π(·). Unlike classical Monte

Carlo estimators, Jˆm is not based on iid random elements. Indeed, the elements of the chain are typically
neither identically distributed nor independent. Given varπ f, the variance of f (·) under the stationary distribution, the accuracy of Jˆm is primarily determined by two factors: (i) the convergence rate of the Markov
chain, and (ii) the correlation between the f (Φk )s when the chain is stationary. These two factors are related,
and can be analyzed jointly under an operator theoretic framework.
The starting point of the operator theoretic approach is the Hilbert space of functions that are square
integrable with respect to the target pdf, π(·). The Markov transition function that gives rise to Φ = {Φk }∞
k=0
defines a linear (Markov) operator on this Hilbert space. (Formal definitions are given in Section 2.) If Φ
is reversible, then it is geometrically ergodic if and only if the corresponding Markov operator admits a
positive spectral gap (Roberts and Rosenthal, 1997; Kontoyiannis and Meyn, 2012). The gap, which is a
real number in (0, 1], plays a fundamental role in determining the mixing properties of the Markov chain,
with larger values corresponding to better performance. For instance, suppose Φ0 has pdf π0 (·) such that
dπ0 /dπ is in the Hilbert space, and let d(Φk ; π) denote the total variation distance between the distribution
of Φk and the chain’s stationary distribution. Then, if δ denotes the spectral gap, we have
d(Φk ; π) ≤ C(1 − δ)k
for all positive integers k, where C depends on π0 but not on k (Roberts and Rosenthal, 1997). Furthermore,
(1 − δ)k gives the maximal absolute correlation between Φj and Φj+k as j → ∞. It follows (see e.g.
√
Mira and Geyer, 1999) that the asymptotic variance of m(Jˆm − J) as m → ∞ is bounded above by
2−δ
varπ f .
δ
Unfortunately, it is impossible to calculate the spectral gaps of the Markov operators associated with practically relevant MCMC algorithms, and even accurately approximating these quantities has proven extremely
difficult. In this paper, we develop a method of estimating the spectral gaps of Markov operators corresponding to a certain class of data augmentation (DA) algorithms (Tanner and Wong, 1987), and then show that
the method can be extended to handle a much larger class of reversible MCMC algorithms.
DA Markov operators are necessarily non-negative. Moreover, any non-negative Markov operator that
is compact has a pure eigenvalue spectrum that is contained in the set [0, 1], and 1 − δ is precisely the
2

second largest eigenvalue. We propose a classical Monte Carlo estimator of 1 − δ for DA Markov operators
that are trace-class, i.e. compact with summable eigenvalues. While compact operators were once thought
to be rare in MCMC problems with uncountable state spaces (Chan and Geyer, 1994), a string of recent
results suggests that trace-class DA Markov operators are not at all rare (see e.g. Qin and Hobert, 2016;
Chakraborty and Khare, 2017; Choi and Román, 2017; Pal et al., 2017). Furthermore, by exploiting a simple
trick, we are able to broaden the applicability of our method well beyond DA algorithms. Indeed, if a
reversible Monte Carlo Markov chain has a Markov transition density (Mtd), and the corresponding Markov
operator is Hilbert-Schmidt, then our method can be utilized to estimate its spectral gap. This is because
the square of such a Markov operator can be represented as a trace-class DA Markov operator. A detailed
explanation is provided in Section 4.
Of course, there is a large literature devoted to developing theoretical bounds on the second largest eigenvalue of a Markov operator (see e.g. Lawler and Sokal, 1988; Sinclair and Jerrum, 1989; Diaconis and Stroock,
1991). However, these results are typically not useful in situations where the state space, S, is uncountable
and multi-dimensional, which is our main focus. There also exist computational methods for approximating the eigenvalues of a Hilbert-Schmidt operator (see e.g. Koltchinskii and Giné, 2000; Ahues et al., 2001,
§4.2). Unfortunately, these methods require a closed form kernel function, which is typically not available
in the MCMC context. There are still other methods based on simulation. Most notably, Garren and Smith
(2000) used simulations of a reversible chain to estimate the second largest eigenvalue of its operator (assuming it’s Hilbert-Schmidt). Their approach is reminiscent of the so-called power method from computer
science, and we use these ideas as well. The main difference between their method and ours is that we exploit
the specific structure of the Mtd associated with the DA algorithm. This makes our method much simpler to
implement computationally, and our results easier to interpret. The power of our method comes at the price
of being computationally intensive, especially when the target posterior is based on a large sample.
The rest of the paper is organized as follows. The notion of Markov operator is formalized in Section 2.
In Section 3, it is shown that the second largest eigenvalue of a non-negative trace-class operator can be
bounded above and below by functions of the power sums of the operator’s eigenvalues. In Section 4, DA
Markov operators are formally defined, and the sum of the kth (k ∈ N) power of the eigenvalues of a
trace-class DA Markov operator is related to a functional of its Mtd. This functional is usually a multidimensional integral, and a classical Monte Carlo estimator of it is developed in Section 5. The efficiency
of the Monte Carlo estimator is studied in Section 6. Finally, in Section 7 we apply our methods to a few
well-known MCMC algorithms. Our examples include Albert and Chib’s (1993) DA algorithm for Bayesian
probit regression, and a DA algorithm for Bayesian linear regression with non-Gaussian errors (Liu, 1996).

3

2 Markov operators
Assume that the Markov chain Φ has a Markov transition density, p(u, ·), u ∈ S, such that, for any measurable A ⊂ S and u ∈ S,
P(Φk ∈ A|Φ0 = u) =
where

Z

p(k) (u, u′ ) µ(du′ ) ,
A


p(u, ·)
k=1
p(k) (u, ·) := R
 p(k−1) (u, u′ )p(u′ , ·) µ(du′ ) k > 1
S

is the k-step Mtd corresponding to p(u, ·). We will assume throughout that Φ is Harris ergodic, i.e. irreducible, aperiodic and Harris recurrent. Define a Hilbert space consisting of complex valued functions on S
that are square integrable with respect to π(·), namely
Z
o

|f (u)|2 π(u) µ(du) < ∞ .
L (π) := f : S → C 
2

n

S

For f, g ∈ L2 (π), their inner product is given by
hf, giπ =

Z

f (u)g(u)π(u) µ(du) .
S

We assume that U is countably generated, which implies that L2 (π) is separable and admits a countable
orthonormal basis (see e.g. Billingsley, 1995, Theorem 19.2). The transition density p(u, ·), u ∈ S defines

the following linear operator P. For any f ∈ L2 (π),
P f (u) =

Z

p(u, u′ )f (u′ ) µ(du′ ) .
S

The spectrum of a linear operator L is defined to be

σ(L) = λ ∈ C (L − λI)−1 doesn’t exist or is unbounded ,

where I is the identity operator. It is well-known that σ(P ) is a closed subset of the unit disk in C. Let
f0 ∈ L2 (π) be the normalized constant function, i.e. f0 (u) ≡ 1, then P f0 = f0 . (This is just a fancy way of
saying that 1 is an eigenvalue with constant eigenfunction, which is true of all Markov operators defined by
ergodic chains.) Denote by P0 the operator such that P0 f = P f − hf, f0 iπ f0 for all f ∈ L2 (π). Then the
spectral gap of P is defined as
n
o
δ = 1 − sup |λ| λ ∈ σ(P0 ) .

For the remainder of this section, we assume that P is non-negative (and thus self-adjoint) and compact.
This implies that σ(P ) ⊂ [0, 1], and that any non-vanishing element of σ(P ) is necessarily an eigenvalue
4

of P . Furthermore, there are at most countably many eigenvalues, and they can accumulate only at the
origin. Let λ0 , λ1 , . . . , λκ be the decreasingly ordered strictly positive eigenvalues of P taking into account
multiplicity, where 0 ≤ κ ≤ ∞. Then λ0 = 1 and λ1 is what we previously referred to as the “second largest
eigenvalue” of the Markov operator. If κ = 0, we set λ1 = 0 (which corresponds to the trivial case where
{Φk }∞
k=0 are iid). Since Φ is Harris ergodic, λ1 must be strictly less than 1. Also, the compactness of P
implies that of P0 , and it’s easy to show that σ(P0 ) = σ(P )\{1}. Hence, Φ is geometrically ergodic and the
spectral gap is
δ = 1 − λ1 > 0 .
For further background on the spectrum of a linear operator, see e.g. Helmberg (2014) or Ahues et al. (2001).

3 Power sums of eigenvalues
We now develop some results relating λ1 to the power sum of P ’s eigenvalues. We assume throughout this
section that P is non-negative and trace-class (compact with summable eigenvalues). For any positive integer
k, let
sk =

κ
X

λki ,

i=0

and define s0 to be infinity. The first power sum, s1 , is the trace norm of P (see e.g. Conway, 1990, 2000),
√
while s2 is the Hilbert-Schmidt norm of P. That P is trace-class implies s1 < ∞, and it’s clear that sk is
decreasing in k.
The magnitude of sk is directly related to the convergence behavior of the chain. For instance, suppose
that the chain starts at a point mass Φ0 = u, then the chi-square distance between the distribution of Φk and
the stationary distribution is given by (see e.g. Diaconis et al., 2008)
χ2k (u) :=

Z

SU

2
κ
X
p(k) (u, u′ ) − π(u′ )
2
′
λ2k
µ(du ) =
i |fi (u)| ,
π(u′ )
i=1

where fi : SU → C is the normalized eigenfunction corresponding to λi . It follows that
s2k :=

κ
X
i=1

λ2k
i

=

Z

SU

χ2k (u)π(u) µ(du),

which is the average of χ2k (u) under π. More importantly, one can use functions of sk to bound λ1 , and thus
the spectral gap.
Observe that,
λ1 ≤ uk := (sk − 1)1/k ,
5

∀k ∈ N.

Moreover, if κ ≥ 1, then it’s easy to show that
λ1 ≥ lk :=

sk − 1
,
sk−1 − 1

∀k ∈ N.

We now show that, in fact, these bounds are monotone in k and converge to λ1 .
Proposition 1. As k → ∞,
uk ↓ λ1 ,

(1)

lk ↑ λ1 .

(2)

and if furthermore κ ≥ 1,

Proof. We begin with (1). When κ = 0, sk ≡ 1 and the conclusion follows. Suppose κ ≥ 1, and that the
second largest eigenvalue is of multiplicity m, i.e.
1 = λ0 > λ1 = λ2 = · · · = λm > λm+1 ≥ · · · ≥ λκ > 0.
If κ = m, then sk − 1 = mλk1 for all k ≥ 1 and the proof is trivial. Suppose for the rest of the proof that
P
κ ≥ m + 1. For positive integer k, let rk = κi=m+1 λki < ∞. Then rk > 0, and
rk+1
≤ λm+1 < λ1 .
rk

Hence,

k−1
r1 λm+1
rk
rk
≤ lim
= 0.
lim
= lim
k→∞ mλk
k→∞ sk − 1 − rk
k→∞ mλk
1
1

It follows that
log uk = log λ1 +
Finally,
1/(k+1)

uk+1 < λ1

κ
X
i=1

and (1) follows.

λki

1
1
log m + log(1 + o(1)) → log λ1 .
k
k

1/(k+1)

≤

κ
X
i=1

λki

κ
1/[k(k+1)]  X

Now onto (2). We have already shown that
sk − 1 = mλk1 (1 + o(1)) .
Thus,
lk =

mλk1 (1 + o(1))
→ λ1 .
mλ1k−1 (1 + o(1))

6

i=1

λki

1/(k+1)

= uk ,

To show that lk is increasing in k, which would complete the proof, we only need note that
(sk+1 − 1)(sk−1 − 1) =
=
≥

κ
X

λk+1
i

κ
X

λjk−1

j=1
κ
κ
1 X X k−1 k−1 2
λi λj (λi
2
i=1 j=1
κ
κ X
X
λki λkj
i=1 j=1
i=1

+ λ2j )

= (sk − 1)2 .

Suppose now that we are interested in the convergence behavior of a particular Markov operator that
is known to be non-negative and trace-class. If it is possible to estimate sk , then Proposition 1 provides a
method of getting approximate bounds on λ1 . When a DA Markov operator is trace-class, there is a nice
integral representation of sk that leads to a simple, classical Monte Carlo estimator of sk . In the following
section, we describe some theory for DA Markov operators, and in Section 5, we develop a classical Monte
Carlo estimator of sk .

4 Data augmentation operators and an integral representation of sk
In order to formally define DA, we require a second measure space. Let (SV , V, ν) be a σ-finite measure
space such that V is countably generated. Also, rename S and π, SU and πU , respectively. Consider the
random element (U, V ) taking values in SU × SV with joint pdf πU,V (·, ·). Suppose the marginal pdf of U is
the target, πU (·), and denote the marginal pdf of V by πV (·). We further assume that the conditional densities
πU |V (u|v) := πU,V (u, v)/πV (v) and πV |U (v|u) := πU,V (u, v)/πU (u) are well defined almost everywhere
in SU × SV . Recall that Φ is a Markov chain on the state space SU with Mtd p(u, ·), u ∈ SU . We call Φ a
DA chain, and accordingly, P a DA operator, if p(u, ·) can be expressed as
p(u, ·) =

Z

SV

πU |V (·|v)πV |U (v|u) ν(dv) .

(3)

Such a chain is necessarily reversible with respect to πU (·). To simulate it, in each iteration, one first draws
the latent element V using πV |U (·|u), where u ∈ SU is the current state, and then given V = v, one updates
the current state according to πU |V (·|v). A DA operator is non-negative, and thus possesses a positive
spectrum (Liu et al., 1994).

7

Assume that (3) holds. Given k ∈ N, the power sum of P ’s eigenvalues, sk , if well defined, is closely

related to the diagonal components of p(k) (·, ·). Just as we can calculate the sum of the eigenvalues of a
R
matrix by summing its diagonals, we can obtain sk by evaluating SU p(k) (u, u) µ(du). Here is a formal
statement.

Theorem 1. The DA operator P is trace-class if and only if
Z
p(u, u) µ(du) < ∞ .

(4)

SU

If (4) holds, then for any positive integer k,
sk :=

κ
X
i=0

λki

=

Z

p(k) (u, u) µ(du) .

(5)

SU

Theorem 1 is a combination of a few standard results in classical functional analysis. It is fairly wellknown, but we were unable to find a complete proof in the literature. An elementary proof is given in the
appendix for completeness. For a more modern version of the theorem, see Brislawn (1988).
It is often possible to exploit Theorem 1 even when Φ is not a DA Markov chain. Indeed, suppose
that Φ is reversible, but is not a DA chain. Then P is not a DA operator, but P 2 , in fact, is. (Just take
πU,V (u, v) = πU (u)p(u, v).) If, in addition, P is Hilbert-Schmidt, which is equivalent to
Z Z
(p(u, u′ ))2 πU (u)
µ(du) µ(du′ ) < ∞ ,
πU (u′ )
SU SU
then by a simple spectral decomposition (see e.g. Helmberg, 2014, §28 Corollary 2.1) one can show that P 2
is trace-class, and its eigenvalues are precisely the squares of the eigenvalues of P . In this case, the spectral
gap of P can be expressed as 1 minus the square root of P 2 ’s second largest eigenvalue. Moreover, by
R
Theorem 1, for k ∈ N, the sum of the kth power of P 2 ’s eigenvalues is equal to SU p(2k) (u, u) µ(du) < ∞.
We now briefly describe the so-called sandwich algorithm, which is a variant of DA that involves an

extra step sandwiched between the two conditional draws of DA (Liu and Wu, 1999; Hobert and Marchev,
2008). Let s(v, ·), v ∈ SV be a Markov transition function (Mtf) with invariant density πV (·). Then
Z Z
πU |V (·|v ′ )s(v, dv ′ )πV |U (v|u)ν(dv)
p̃(u, ·) =
SV

(6)

SV

is an Mtd with invariant density πU (·). This Mtd defines a new Markov chain, call it Φ̃, which we refer
to as a sandwich version of the original DA chain, Φ. To simulate Φ̃, in each iteration, the latent element is first drawn from πV |U (·|u), and then updated using s(v, ·) before the current state is updated ac-

cording to πU |V (·|v ′ ). Sandwich chains often converge much faster than their parent DA chains (see e.g.
Khare and Hobert, 2011).
8

Of course, p̃(u, ·) defines a Markov operator on L2 (πU ), which we refer to as P̃ . It is easy to see that,
if the Markov chain corresponding to s(v, ·) is reversible with respect to πV (·), then p̃(u, ·) is reversible

with respect to πU (·). Thus, when s(v, ·) is reversible, P̃ 2 is a DA operator. Interestingly, it turns out that
p̃(u, ·) can often be re-expressed as the Mtd of a DA chain, in which case P̃ itself is a DA operator. Indeed,
a sandwich Mtd p̃(u, ·) is said to be “representable” if there exists a random element Ṽ in SV such that
′

p̃(u, u ) =

Z

SV

πU |Ṽ (u′ |v)πṼ |U (v|u) ν(dv) ,

(7)

where πU |Ṽ (u′ |v) and πṼ |U (v|u) have the apparent meanings (see, e.g. Hobert, 2011). It is shown in Proposition 2 in Section 5 that when P is trace-class and p̃(u, ·) is representable, P̃ is also trace-class. In this case,

let {λ̃i }κ̃i=0 be the decreasingly ordered positive eigenvalues of P̃ taking into account multiplicity, where
0 ≤ κ̃ ≤ ∞. Then λ̃0 = 1, and λ̃1 ≤ λ1 < 1 (Hobert and Marchev, 2008). For a positive integer k, we will
P
denote κ̃i=0 λ̃ki by s̃k . Henceforth, we assume that p̃(u, ·) is representable and we treat P̃ as a DA operator.
R
It follows from Theorem 1 that in order to find sk or s̃k , all we need to do is evaluate SU p(k) (u, u) µ(du)
R
or SU p̃(k) (u, u) µ(du), where p̃(k) (u, ·) is the k-step Mtd of the sandwich chain. Of course, calculating

these integrals (in non-toy problems) is nearly always impossible, even for k = 1. In the next section, we
introduce a method of estimating these two integrals using classical Monte Carlo.
Throughout the remainder of the paper, we assume that P is a DA operator with Mtd given by (3), and
that (4) holds.

5 Classical Monte Carlo
Consider the Mtd given by
a(u, ·) =

Z

SV

Z

SV

πU |V (·|v ′ )r(v, dv ′ )πV |U (v|u) ν(dv) ,

(8)

where r(v, ·), v ∈ SV is an Mtf on SV with invariant pdf πV (·). We will show in this section that this
form has utility beyond constructing sandwich algorithms. Indeed, the k-step Mtd of a DA algorithm (or
a sandwich algorithm) can be re-expressed in the form (8). This motivates the development of a general
R
method for estimating the integral SU a(u, u) µ(du), which is the main topic of this section.

We begin by showing how p(k) (u, ·), u ∈ SU can be written in the form (8). The case k = 1 is trivial.

Indeed, if r(v, ·) is taken to be the kernel of the identity operator, then a(u, ·) = p(u, ·). Define an Mtd
q(v, ·), v ∈ SV by
q(v, ·) =

Z

SU

πV |U (·|u)πU |V (u|v) µ(du),

9

and let q (k) (v, ·), k ≥ 1 denote the corresponding k-step Mtd. If we let
r(v, dv ′ ) = q (k−1) (v, v ′ ) ν(dv ′ ), v ∈ SV
for k ≥ 2, then a(u, ·) = p(k) (u, ·). Next, consider the sandwich Mtd p̃(k) (u, ·), u ∈ SU . Again, the k = 1
case is easy. Taking
r(v, ·) = s(v, ·)
yields a(u, ·) = p̃(u, ·). Now let
q̃(v, ·) =

Z

SU

Z

SV

s(v ′ , ·)πV |U (v ′ |u)πU |V (u|v) ν(dv ′ ) µ(du) ,

and denote the corresponding k-step transition function by q̃ (k) (v, ·). Then taking
r(v, ·) =

Z

SV

q̃ (k−1) (v ′ , ·)s(v, dv ′ )

when k ≥ 2 yields a(u, ·) = p̃(k) (u, ·).
The following proposition shows that, when P is trace-class,
R

Proposition 2.
Proof. That
Z

R

SU

SU

a(u, u) µ(du) < ∞.

SU

R

SU

a(u, u) µ(du) is finite.

a(u, u) µ(du) < ∞ is equivalent to


πU,V (u, v)
πU,V (u, v ′ )
′
r(v, dv )
πU (u)πV (v) ν(dv) µ(du) < ∞.
πU (u)πV (v ′ )
πU (u)πV (v)

Z

Z

SV

SV

Note that
Z

SU



πU,V (u, v)
πU (u)πV (v)

2

πU (u)πV (v) µ(du)ν(dv) =

Z

SU

p(u, u) µ(du) < ∞.

(9)

(10)

and by Jensen’s inequality,
Z

SU

SV

≤

Z

Z

=

Z

Z

=

SU

ZSU

2
πU,V (u, v ′ )
′
r(v, dv ) πU (u)πV (v) ν(dv) µ(du)
′
SV πU (u)πV (v )

Z 
πU,V (u, v ′ ) 2
r(v, dv ′ )πU (u)πV (v) ν(dv) µ(du)
′)
π
(u)π
(v
U
V
SV


′
πU,V (u, v ) 2
πU (u)πV (v ′ ) ν(dv ′ ) µ(du)
πU (u)πV (v ′ )

Z

Z

SV

SV

p(u, u) µ(du)

SU

< ∞.
The inequality (9) follows from (10), (11), and the Cauchy-Schwarz inequality.
10

(11)

Combining Proposition 2 and Theorem 1 leads to the following result: If P is trace-class and p̃(u, ·)
is representable, then P̃ is also trace-class. This is a generalization of Khare and Hobert’s (2011) Theorem 1, which states that, under a condition on s(v, dv ′ ) that implies representability, the trace-class-ness of
P implies that of P̃ .
We now develop a classical Monte Carlo estimator of

R

SU

a(u, u) µ(du). Let ω : SV → [0, ∞) be a pdf

that is almost everywhere positive. We will exploit the following representation of the integral of interest:
 Z
Z Z 
Z

πV |U (v|u) 
(12)
πU |V (u|v ′ )r(v, dv ′ ) ω(v) µ(du) ν(dv) .
a(u, u) µ(du) =
ω(v)
SV
SV SU
SU
Clearly,
η(u, v) :=
(U ∗ , V ∗ )

defines a pdf on SU × SV , and if
Z

SU

Z

SV


πU |V (u|v ′ )r(v, dv ′ ) ω(v)

has joint pdf η(·, ·), then


πV |U (V ∗ |U ∗ )
a(u, u) µ(du) = E
.
ω(V ∗ )

Therefore, if {(Ui∗ , Vi∗ )}N
i=1 are iid random elements from η(·, ·), then
N
1 X πV |U (Vi∗ |Ui∗ )
N
ω(Vi∗ )

(13)

i=1

is a strongly consistent and unbiased estimator of
central to our discussion.

R

SU

a(u, u) µ(du). This is the Monte Carlo formula that is

Of course, we are mainly interested in the cases a(u, ·) = p(k) (u, ·) or a(u, ·) = p̃(k) (u, ·). We now

develop algorithms for drawing from η(·, ·) in these two situations. First, assume a(u, ·) = p(k) (u, ·). If
k = 1, then r(u, ·) is the kernel of the identity operator, and
η(u, v) = πU |V (u|v)ω(v) .
If k ≥ 2, then r(v, dv ′ ) = q (k−1) (v, v ′ ) dv ′ , and

Z
Z
′ (k−1)
′
′
πU |V (u|v )q
(v, v ) ν(dv ) ω(v) =
η(u, v) =

SU

SV


p(k−1) (u′ , u)πU |V (u′ |v) µ(du′ ) ω(v) .

Thus, when k ≥ 2, we can draw from η(u, v) as follows: Draw V ∗ ∼ ω(·), then draw U ′ ∼ πU |V (·|v ∗ ), then
draw U ∗ ∼ p(k−1) (u′ , ·), and return (u∗ , v ∗ ). Of course, we can draw from p(k−1) (u′ , ·) by simply running

k − 1 iterations of the original DA algorithm from starting value u′ . We formalize all of this in Algorithm 1.

Algorithm 1: Drawing (U ∗ , V ∗ ) ∼ η(u, v) µ(du) ν(dv) when a(u, ·) = p(k) (u, ·).
11

1. Draw V ∗ from ω(·).
2. Given V ∗ = v ∗ , draw U ′ from πU |V (·|v ∗ ).
3. If k = 1, set U ∗ = U ′ . If k ≥ 2, given U ′ = u′ , draw U ∗ from p(k−1) (u′ , ·) by running k − 1 iterations
of the DA algorithm.

Similar arguments lead to the following algorithm for the sandwich algorithm

Algorithm 1S: Drawing (U ∗ , V ∗ ) ∼ η(u, v) µ(du) ν(dv) when a(u, ·) = p̃(k) (u, ·)
1. Draw V ∗ from ω(·).
2. Given V ∗ = v ∗ , draw V ′ from s(v ∗ , ·).
3. Given V ′ = v ′ draw U ′ from πU |V (·|v ′ ).
4. If k = 1, set U ∗ = U ′ . If k ≥ 2, given U ′ = u′ , draw U ∗ from p̃(k−1) (u′ , ·) by running k − 1 iterations
of the sandwich algorithm.
It is important to note that we do not need to know the representing conditionals πU |Ṽ (·|v) and πṼ |U (·|u)
from (7) in order to run Algorithm 1S.
As with all classical Monte Carlo techniques, the key to successful implementation is a finite variance.
Define

πV |U (V ∗ |U ∗ )
D = var
ω(V ∗ )
2

Of course, D 2 < ∞ if and only if
Z Z
SV

SU





πV |U (v|u)
ω(v)

2



.

η(u, v) µ(du) ν(dv) < ∞ .

(14)

The following theorem provides a sufficient condition for finite variance.
Theorem 2. The variance, D 2 , is finite if
Z

SV

Z

SU

πV3 |U (v|u)πU |V (u|v)
ω 2 (v)

µ(du) ν(dv) < ∞.

Proof. First, note that (14) is equivalent to
Z Z  π 2 (v|u)  R π (u|v ′ )r(v, dv ′ ) 
V |U
SV U |V
SV

SU

πV (v)ω(v)

πU (u)

12

πU (u)πV (v) µ(du) ν(dv) < ∞.

(15)

Now, it follows from (15) that
Z Z
SV

 π2

V |U (v|u)

SU

Moreover, by Jensen’s inequality,
Z Z R

SV

SV

Z

Z

=

Z

Z

SV

ZSV

SU

SU

πU (u)πV (v) µ(du) ν(dv) < ∞.

πU |V (u|v ′ )r(v, dv ′ )
πU (u)

SU

≤

=

πV (v)ω(v)

2

2

πU (u)πV (v) µ(du) ν(dv)


πU |V (u|v ′ ) 2
r(v, dv ′ )πU (u)πV (v) µ(du) ν(dv)
π
(u)
U
SV


πU |V (u|v ′ ) 2
πU (u)πV (v ′ ) µ(du) ν(dv ′ )
πU (u)

Z

(16)



(17)

p(u, u) µ(du)

SU

< ∞.
The conclusion now follows from (16), (17), and Cauchy-Schwarz.
Theorem 2 implies that an ω(·) with heavy tails is more likely to result in finite variance (which is
not surprising). It might seem natural to take ω(·) = πV (·). However, in practice, we are never able to
draw from πV (·). (If we could do that, we would not need MCMC.) Moreover, in Section 7 we provide
an example where taking ω(·) to be πV (·) leads to an infinite variance, whereas a heavier-tailed alternative
leads to a finite variance. On the other hand, it can be beneficial to use ω(·)s resembling πV (·), as we argue
in Section 6.
Let ψ : SU → [0, ∞) be a pdf that is positive almost everywhere. The following dual of (12) may also
R
be used to represent SU a(u, u) µ(du):
Z Z Z
Z
πU |V (u|v)
r(v ′ , dv)πV |U (v ′ |u)ψ(u) ν(dv ′ ) µ(du) .
a(u, u) µ(du) =
ψ(u)
SU SV SV
SU
Now suppose that {(Ui∗ , Vi∗ )}N
i=1 are iid from

Z
′
′
′
r(v , dv)πV |U (v |u) ν(dv ) ψ(u) µ(du) .
ζ(u, v) µ(du) ν(dv) =
SV

The analogue of (13) is the following classical Monte Carlo estimator of
N
1 X πU |V (Ui∗ |Vi∗ )
.
N
ψ(Ui∗ )

R

SU

a(u, u) µ(du):

i=1

We now state the obvious analogues of Algorithms 1 and 1S.

Algorithm 2: Drawing (U ∗ , V ∗ ) ∼ ζ(u, v) µ(du) ν(dv) when a(u, ·) = p(k) (u, ·).
13

(18)

1. Draw U ∗ from ψ(·).
2. If k = 1, set U ′ = U ∗ . If k ≥ 2, given U ∗ = u∗ , draw U ′ from p(k−1) (u∗ , ·).
3. Given U ′ = u′ , draw V ∗ from πV |U (·|u′ ).

Algorithm 2S: Drawing (U ∗ , V ∗ ) ∼ ζ(u, v) µ(du) ν(dv) when a(u, ·) = p̃(k) (u, ·).
1. Draw U ∗ from ψ(·).
2. If k = 1, set U ′ = U ∗ . If k ≥ 2, given U ∗ = u∗ , draw U ′ from p̃(k−1) (u∗ , ·).
3. Given U ′ = u′ , draw V ′ from πV |U (·|u′ ).
4. Given V ′ = v ′ , draw V ∗ from s(v ′ , ·).
Let D ′2 be the variance of πU |V (U ∗ |V ∗ )/ψ(U ∗ ) under ζ. To ensure that it’s finite, we only need
Z

SU

Z

SV



Z

SV

πU |V (u|v)
ψ(u)

2

r(v ′ , dv)πV |U (v ′ |u)ψ(u) ν(dv ′ ) µ(du) < ∞ .

(19)

The following result is the analogue of Theorem 2.
Corollary 1. The variance, D ′2 , is finite if
Z

SU

Z

πU3 |V (u|v)πV |U (v|u)
SV

ψ 2 (u)

ν(dv) µ(du) < ∞ .

(20)

Proof. Note that the left hand side of (19) is equal to
Z

SU

Z

SV

Z

πU2 |V (u|v)

SV



πV |U (v ′ |u)
r(v ′ , dv)
πU (u)πV (v ′ ) ν(dv ′ ) µ(du) .
ψ(u)πU (u)
πV (v ′ )

Apply the Cauchy-Schwarz inequality, then utilize Jensen’s inequality to get rid of r(v ′ , dv), and finally
make use of (20) and the fact that P is trace-class.
Suppose that we have obtained estimates of sk and sk−1 based on (13) or (18), call them s∗k and s∗k−1 .
Then u∗k = (s∗k − 1)1/k and lk∗ = (s∗k − 1)/(s∗k−1 − 1) serve as point estimates of uk and lk , respectively.
When our estimators have finite variances, we can acquire, via the delta method, confidence intervals for uk
and lk . Assume that a confidence interval for lk is (ak , bk ) and a confidence interval for uk is (ck , dk ), then
(ak , dk ) is an interval estimate for λ1 . Interval estimates of λ̃1 can be derived in a similar fashion.

14

6 Efficiency of the algorithm
To obtain an interval estimate of λ1 based on (13) or (18), one needs to run N iterations of Algorithm 1
or 2. If the time needed to simulate one step of the DA chain is τ , then the time needed to run N iterations
of Algorithm 1 or 2 is approximately kN τ . Given k and N , the accuracy of the estimate depends on two
factors: 1. The distance between lk and uk , and 2. The errors in the estimates, lk∗ and u∗k . We now briefly
analyze these two factors, and give some additional guidelines regarding the choice of ω(·) and ψ(·).
As before, suppose that
1 = λ0 > λ1 = λ2 = · · · = λm > λm+1 ≥ · · · ≥ λκ > 0
for some m < ∞. Clearly, as k → ∞,



sk − 1 = λk1 m + O λkm+1 /λk1 .
Hence, as k → ∞,
lk :=
and




sk − 1
k−1
/λ1k−1 ,
= λ1 1 + O λm+1
sk−1 − 1

uk := (sk − 1)1/k



= λ1 m1/k 1 + O k−1 λkm+1 /λk1


 
1 + O k−1 λkm+1 /λk1 .
= λ1 1 + (log m)O k−1

λ1 1 + O k−1 λk /λk  m = 1
2
1
=

λ 1 + O k−1
m > 1.
1

Depending on whether m = 1 or not, uk − lk decreases at either a geometric or polynomial rate as k grows.
The errors of lk∗ and u∗k arise from those of s∗k and s∗k−1 . We now consider the estimator (18) for estimating

sk . Its variance is given by
(Z Z Z 
)

πU |V (u|v) 2 ′
D ′2
1
′
′
2
=
r(v , dv)πV |U (v |u)ψ(u) ν(dv ) µ(du) − sk
N
N
ψ(u)
SU SV SV
(Z Z Z
)
πU2 |V (u|v)
1
=
πV |U (v|u′ )p(k) (u, u′ ) µ(du′ ) ν(dv) µ(du) − s2k .
N
ψ(u)
SU SV SU
Note that

pk (u, v), (u′ , v ′ ) := πV |U (v ′ |u′ )p(k) (u, u′ )

15

gives the k-step Mtd of a Gibbs chain whose stationary pdf is πU,V (·, ·). Thus, under suitable conditions, for
almost any u ∈ SU ,
lim sk (u) := lim

k→∞

Z

k→∞ SV

=

Z

Z

Z

πU2 |V (u|v)
SU

πU2 |V (u|v)

ψ(u)
p(u, u)πU (u)
.
=
ψ(u)
SV

πV |U (v|u′ )p(k) (u, u′ ) µ(du′ ) ν(dv)

ψ(u)

πU,V (u′ , v) µ(du′ ) ν(dv)

SU

As k → ∞, we expect
′2

D =

Z

SU

sk (u) µ(du) −

s2k

→

Z

SU

p(u, u)πU (u)
µ(du) − 1.
ψ(u)

Suppose that ψ(u) ≈ πU (u), then heuristically,
Z
Z
p(u, u)πU (u)
p(u, u)µ(du) − 1 = s1 − 1.
µ(du) − 1 ≈
ψ(u)
SU
SU
Thus, if the sum of P ’s eigenvalues, s1 , is relatively small, we recommend picking ψ(·)s that resemble πU (·),
with possibly heavier tails (to ensure that the moment condition (20) holds). By a similar argument, when
using the estimator (13), picking ω(·)s that resemble πV (·) is likely to control D 2 around s1 − 1 for large ks.
While (under suitable conditions) the variance of s∗k converges to a constant as k → ∞, this is not the

case for u∗k and lk∗ (because uk and lk are non-linear in sk and sk−1 ). In fact, using the delta method, one
can show that these variances are unbounded. Thus, there’s a trade-off between decreasing uk − lk (by

increasing k) and controlling the errors of u∗k and lk∗ . We do not recommend increasing k indefinitely. As

long as k is large enough so that sk − 1 is significantly smaller than 1, uk serves as a non-trivial (and often
decent) upper bound for λ1 .

7 Examples
In this section, we apply our Monte Carlo technique to several common Markov operators. In particular, we
examine one toy Markov chain, and two practically relevant Monte Carlo Markov chains. In the two real
examples, we are able to take advantage of existing trace-class proofs to establish that (15) (or (20)) hold for
suitable ω(·) (or ψ(·)).

7.1 Gaussian chain
We begin with a toy example. Let SU = SV = R, πU (u) ∝ exp(−u2 ), and
n

u 2 o
πV |U (v|u) ∝ exp − 4 v −
.
2
16

Then
πU |V (u|v) ∝ exp{−2(u − v)2 } .
This leads to one of the simplest DA chains known. Indeed, the Mtd,
p(u, ·) =

Z

πU |V (·|v)πV |U (v|u) dv ,

R

can be evaluated in closed form, and turns out to be a normal pdf. The spectrum of the corresponding Markov
operator, P , has been studied thoroughly (see e.g. Diaconis et al., 2008). It is easy to verify that (4) holds,
so P is trace-class. In fact, κ = ∞, and for any non-negative integer i, λi = 1/2i . Thus, the second largest
eigenvalue, λ1 , and the spectral gap, δ, are both equal to 1/2. Moreover, for any positive integer k,
∞
X
1
1
=
.
sk =
ik
2
1 − 2−k
i=0

We now pretend to be unaware of this spectral information, and use (13) to estimate {sk , lk , uk }4k=1 .
Recall that lk and uk are lower and upper bounds for λ1 , respectively. Note that
Z

R

πV3 |U (v|u)πU |V

6 2
(u|v) du ∝ exp − v .
5


It follows that, if we take ω(v) ∝ exp(−v 2 /2), then (15) holds, and our estimator of sk has finite variance.

We used a Monte Carlo sample size of N = 1 × 105 to form our estimates, and the results are shown in
Table 1.
Table 1: Estimated power sums of eigenvalues for the Gaussian chain
k

Est. sk

√
Est. D/ N

Est. lk

Est. uk

1

1.996

0.004

0.000

0.996

2

1.331

0.004

0.333

0.575

3

1.142

0.004

0.429

0.522

4

1.068

0.004

0.482

0.511

Note that the estimates of the sk s are quite good. We constructed 95% confidence intervals (CIs) for l4
and u4 via the delta method, and the results were (0.442, 0.522) and (0.498, 0.524), respectively.
Remark 1. It is interesting that, if ω(·) is set to be πV (·), which seems natural, then (15) fails to hold. In
fact, the estimator of s1 has infinite variance in this case. Indeed, recall that the estimator given in (13) has

17

the form
N
1 X πV |U (Vi∗ |Ui∗ )
,
N
ω(Vi∗ )
i=1

where

(Ui∗ , Vi∗ )

are iid, and

(U1∗ , V1∗ )

has pdf given by

η(u, v) = πU |V (u|v)ω(v).
Hence, the variance of the estimator, D 2 , is finite if and only if
Z

SU

Z

πV2 |U (v|u)πU |V (u|v)
ω(v)

SV

dv du < ∞.

(21)

If ω(·) = πV (·) ∝ exp(−2v 2 ), then the left hand side of (21) becomes proportional to
Z Z
Z Z
o
n

u 2
exp(−8v 2 + 12uv − 4u2 ) dv du,
− 2(u − v)2 + 2v 2 dv du =
exp − 8 v −
2
SU SV
SU SV
which is infinite.

7.2 Bayesian probit regression
Let Y1 , Y2 , . . . , Yn be independent Bernoulli random variables with P(Y1 = 1|β) = Φ(xTi β), where xi , β ∈
Rp , and Φ(·) is the cumulative distribution function of the standard normal distribution. Take the prior on

β to be Np (Q−1 w, Q−1 ), where w ∈ Rp and Q is positive definite. The resulting posterior distribution is

intractable, but Albert and Chib (1993) devised a DA algorithm to sample from it. Let z = (z1 , z2 , . . . , zn )T

be a vector of latent variables, and let X be the design matrix whose ith row is xTi . The Mtd of the Albert
and Chib (AC) chain, p(β, ·), β ∈ Rp , is characterized by


−1
 T

−1

1
T
T
T
T
T
w +X z
X X +Q β − X X +Q
w +X z
,
πU |V (β|z) ∝ exp − β − X X +Q
2
and

πV |U (z|β) ∝

n
Y
i=1

exp




2

1
T
− zi − xi β
IR+ (yi − 0.5)zi .
2

The first conditional density, πU |V (·|z), is a multivariate normal density, and the second conditional density,
πV |U (·|β), is a product of univariate truncated normal pdfs.
A sandwich step can be added to facilitate the convergence of the AC chain. Chakraborty and Khare
(2017) constructed a Haar PX-DA variant of the chain, which is a sandwich chain with transition density of
the form (6) (see also Roy and Hobert (2007)). The sandwich step s(z, dz ′ ) is equivalent to the following
update: z 7→ z ′ = gz, where the scalar g is drawn from the following density:


1 T
T
−1 T
2
T
T
−1
n−1
πG (g|z) ∝ g
exp − z In − X(X X + Q) X zg + z X(X X + Q) wg .
2
18

Note that this pdf is particularly easy to sample from when w = 0.
Chakraborty and Khare showed that P is trace-class when one uses a concentrated prior (corresponding
to Q having large eigenvalues). In fact, the following is shown to hold in their proof.
Proposition 3. Suppose that X is full rank. If all the eigenvalues of Q−1/2 X T XQ−1/2 are less than 7/2,
then for any polynomial function t : Rp → R,
Z
|t(β)|p(β, β) dβ < ∞ .
Rp

We will use the estimator (18). The following proposition provides a class of ψ(·)s that lead to estimators
with finite variance.
Proposition 4. Suppose the hypothesis in Proposition 3 holds. If ψ(·) is the pdf of a p-variate t-distribution,
i.e.
ψ(β) ∝

−(a+p)/2

1
1 + (β − b)T Σ−1 (β − b)
a

for some b ∈ Rp , positive definite matrix Σ ∈ Rp×p , and positive integer a, then the estimator (18) has finite
variance.
Proof. Note that for every β and z
πU3 |V (β|z) ≤ CπU |V (β|z) ,
where C is a constant. Hence, by Proposition 3, for any polynomial function t : Rp → R,
Z
Z Z
3
|t(β)|p(β, β) dβ < ∞.
|t(β)|πU |V (β|z)πV |U (z|β) dz dβ ≤ C
Rp

Rp

Rn

Since ψ −2 (·) is a polynomial function on Rp , the moment condition (20) holds. The result follows from
Corollary 1.
As a numerical illustration, we apply our method to the Markov operator associated with the AC chain
corresponding to the famous “lupus data” of van Dyk and Meng (2001). In this data set, n = 55 and p = 3.
As in Chakraborty and Khare (2017), we will let w = 0 and Q = X T X/g, where g = 3.499999. It can be
easily shown that the assumptions in Proposition 3 are met. Chakraborty and Khare compared the AC chain,
Φ, and its Haar PX-DA variant, Φ̃, defined a few paragraphs ago. This comparison was done using estimated
autocorrelations. Their results suggest that Φ̃ outperforms Φ when estimating a certain test function. We go
further and estimate the second largest eigenvalue of each operator.
It can be shown that the posterior pdf, πU (·), is log-concave, and thus possess a unique mode. Let β̂ be the
posterior mode, and Σ̂ the estimated variance of the MLE. We pick ψ(·) to be the pdf of t30 (β̂, (Σ̂−1 +Q)−1 ).
19

This is to say, for any β ∈ Rp ,
o−(p+30)/2
n
1
.
ψ(β) ∝ 1 + (β − β̂)T (Σ̂−1 + Q)(β − β̂)
30
By Proposition 4, this choice of ψ(·) guarantees finite variance. When n is large, ψ(·) is expected to resemble
πU (·). The performance of our method seems insensitive to the degrees of freedom of the t-distribution
(which is set at 30 for illustration).
We used a Monte Carlo sample size of N = 4 × 105 to form our estimates for the DA operator, and
the results are shown in Table 2. Asymptotic 95% CIs for l5 and u5 are (0.397, 0.545) and (0.573, 0.595),
respectively. Using a Bonferroni argument, we can state that asymptotically, with at least 95% confidence,
λ1 ∈ (0.397, 0.595).
Table 2: Estimated power sums of eigenvalues for the AC chain
k

Est. sk

√
Est. D ′ / N

Est. lk

Est. uk

1

6.744

0.072

0.000

5.744

2

2.041

0.007

0.181

1.020

3

1.363

0.004

0.349

0.713

4

1.156

0.004

0.430

0.628

5

1.068

0.003

0.436

0.584

Table 3: Estimated power sums of eigenvalues for the Haar PX-DA version of the AC chain
k

Est. s̃k

√
Est. D ′ / N

Est. ˜lk

Est. ũk

1

3.796

0.012

0.000

1.796

2

1.538

0.004

0.193

0.734

3

1.172

0.004

0.319

0.556

4

1.060

0.003

0.352

0.496

5

1.025

0.003

0.419

0.479

We now consider the sandwich chain, Φ̃. One can show that the Mtd of any Haar PX-DA chain is
representable (Hobert and Marchev, 2008). Hence, P̃ is indeed a DA operator. Recall that {λ̃i }κ̃i=0 , 0 ≤
κ̃ ≤ ∞, denote the decreasingly ordered positive eigenvalues of P̃ . It was shown in Khare and Hobert
20

(2011) that λ̃i ≤ λi for i ∈ N with at least one strict inequality. For a positive integer k,

Pκ̃

k
i=0 λ̃i

is denoted

by s̃k . Let ũk and ˜lk , be the respective counterparts of uk and lk . Estimates of s̃k , k = 1, 2, · · · , 5 using

4 × 105 Monte Carlo samples are given in Table 3. Our estimate of s̃1 − 1 is less than half of s1 − 1,
implying that, in an average sense, the sandwich version of the AC chain reduces the nontrivial eigenvalues
of P by more than half. Asymptotic 95% CIs for l5 and u5 are (0.321, 0.518) and (0.456, 0.503). Thus,
asymptotically, with at least 95% confidence, λ̃1 ∈ (0.321, 0.503). The method does not detect a significant
difference between λ1 and λ̃1 .

7.3 Bayesian linear regression model with non-Gaussian errors
Let Y1 , Y2 , . . . , Yn be independent d-dimensional random vectors from the linear regression model
Yi = β T xi + Σ1/2 εi ,
where xi ∈ Rp is known, while β ∈ Rp×d and the d × d positive definite matrix Σ are to be estimated. The
iid errors, ε1 , ε2 , . . . , εn , are assumed to have a pdf that is a scaled mixture of Gaussian densities:
Z
 u

ud/2
T
exp
−
ε
ε
h(u) du,
fh (ε) =
d/2
2
R+ (2π)
where h(·) is a pdf with positive support, and R+ := (0, ∞). For instance, if d = 1 and h(u) ∝ u−2 e−1/(8u) ,
then ε1 has pdf proportional to e−|ε|/2 .

To perform a Bayesian analysis, we require a prior on the unknown parameter, (β, Σ). We adopt the
(improper) Jeffreys prior, given by 1/|Σ|(d+1)/2 . Let y represent the n × d matrix whose ith row is the
observed value of Yi . The following four conditions, which are sufficient for the resulting posterior to be
proper (Qin and Hobert, 2016; Fernandez and Steel, 1999), will be assumed to hold:
1. n ≥ p + d,
2. (X : y) is full rank, where X is the n × p matrix whose ith row is xTi ,
3.
4.

R

R

R+

ud/2 h(u) du < ∞, and

R+

u−(n−p−d)/2 h(u) du < ∞.

The posterior density is highly intractable, but there is a well-known DA algorithm to sample from it (Liu,
1996). Under our framework, the DA chain Φ is characterized by the Mtd
Z

πU |V (·, ·|z)πV |U (z|β, Σ) dz,
p (β, Σ), (·, ·) =
Rn
+

21

where z = (z1 , z2 , . . . , zn )T ,
πU |V (β, Σ|z) ∝ |Σ|−(n+d+1)/2
πV |U (z|β, Σ) ∝

n
Y

d/2

zi

exp

i=1

n

n
Y
i=1

−

exp

n

−

o
T
zi
yi − β T xi Σ−1 yi − β T xi , and
2

T
o
zi
yi − β T xi Σ−1 yi − β T xi h(zi ) .
2

The first conditional density, πU |V (·, ·|z), characterizes a multivariate normal distribution on top of an inverse
Wishart distribution, i.e. β|Σ, z is multivariate normal, and Σ|z is inverse Wishart. The second conditional
density, πV |U (·|β, Σ), is a product of n univariate densities. Moreover, when h(·) is a standard pdf on R+ ,
these univariate densities are often members of a standard parametric family. The following proposition
about the resulting DA operator is proved in Qin and Hobert (2016).
Proposition 5. Suppose h(·) is strictly positive in a neighborhood of the origin. If there exists ξ ∈ (1, 2)
and δ > 0 such that
Z

0

then P is trace-class.

δ

ud/2 h(u)
du < ∞,
R ξu
d/2 h(v) dv
0 v

When P is trace-class, we can pick an ω(·) and try to make use of (13). A sufficient condition for
the estimator’s variance, D 2 , to be finite is stated in the following proposition, whose proof is given in the
appendix.
Proposition 6. Suppose that h(·) is strictly positive in a neighborhood of the origin. If ω(z) can be written
Q
as ni=1 ωi (zi ), and there exists ξ ∈ (1, 4/3) such that for all i ∈ {1, 2, . . . , n},
Z

R+

u3d/2 h3 (u)
du < ∞,
R ξu
( 0 v d/2 h(v) dv)3 ωi2 (u)

(22)

then (15) holds, and thus by Theorem 2, the estimator (13) has finite variance.
For illustration, take d = 1 and h(u) ∝ u−2 e−1/(8u) . Then ε1 follows a scaled Laplace distribution, and
the model can be viewed as a median regression model with variance Σ unknown. It’s easy to show that
h(·) satisfies the assumptions in Proposition 5, so the resultant DA operator is trace-class. (This result was
actually first proven by Choi and Hobert (2013).) Now let
ω(z) =

n
Y
i=1

ωi (zi ) ∝

n
Y

−3/2 −1/(32zi )

zi

e

.

i=1

The following result shows that this will lead to an estimator with finite variance.

22

Corollary 2. Suppose d = 1, h(u) ∝ u−2 e−1/(8u) , and
ω(z) =

n
Y
i=1

ωi (zi ) ∝

n
Y

zi−α−1 e−γ/zi ,

i=1

where 0 < α < 3/4 and 0 < γ < 3/64. Then the variance, D 2 , is finite.
Proof. In light of Proposition 6, we only need to show that (22) holds for some ξ ∈ (1, 4/3). For any ξ > 0,
making use of the fact that (by monotone convergence theorem)
ξu

Z

lim

u→∞ 0

v 1/2 h(v) dv =

Z

u1/2 h(u) du > 0 ,

R+

one can easily show for any δ > 0,
Z

∞
δ

u3/2 h3 (u)

du =
R ξu
( 0 v 1/2 h(v) dv)3 ωi2 (u)

Z

∞

δ

u2α−5/2 exp{2γ/u − 3/(8u)}
du < ∞.
R ξu
( 0 v 1/2 h(v) dv)3

(23)

On the other hand, using L’Hôpital’s rule, we can see for (1 − 16γ/3)−1 < ξ < 4/3,
lim

u→0

u3/2 h3 (u)
R ξu
( 0 v 1/2 h(v) dv)3 ωi2 (u)

!1/3

u2α/3−5/6 exp{2γ/(3u) − 1/(8u)}
R ξu
−3/2 e−1/(8v) dv
u→0
0 v
n  2γ
1
1 1 o
= lim R(u) exp − −
−
+
u→0
3
8ξ 8 u

= lim

= 0,

where R(u) is a function that is either bounded near the origin or goes to ∞ at the rate of some power
function as u → 0. It follows that for ξ ∈ ((1 − 16γ/3)−1 , 4/3) and small enough δ,
Z

0

δ

u3/2 h3 (u)
du < ∞.
R ξu
( 0 v 1/2 h(v) dv)3 ωi2 (u)

(24)

Combining (23) and (24) yields (22). The result then follows.

We now test the efficiency of the Monte Carlo estimator (13) using some simulated data with d = 1.

23

Here are our simulated X and y:

1
2.32


 1
5.65


 1 −7.69


 1
3.59


 1
5.57

X=
 1 −9.99


 1 −18.88



5.95
 1


 1 −16.39

1
4.75

1 0

2.32



0




1 0 5.65
0



1 0 −7.69
0


0 1
0
3.59 


0 1
0
5.57 


0 1
0
−9.99 


0 1
0
−18.88 



0 0
0
0




0 0
0
0

0 0
0
0



and














y=













0.14
2.99
1.37
−2.55
−3.60
14.86
21.24
2.78
9.14
2.06
















.













The simulation was was based on a linear model containing an intercept, one continuous covariate, a single
factor with three levels, and an interaction between the two. The elements in the second column of X were
independently generated from N (0, 100). Once X was simulated, we generated the data according to
Yi = β∗T xi +

p

Σ ∗ ǫi ,

where β∗ = (1, 0, −0.5, 0.5, 0, −1)T , Σ∗ = 1, {ǫi }10
i=1 are iid, and ǫ1 has pdf given by fh (·) when h(u) ∝

u−2 e−1/(8u) . That is, the errors have a scaled Laplace distribution. Then the DA chain Φ lives in SU =
6
R6 × R+ , and SV is R10
+ . We used a Monte Carlo sample size of N = 2 × 10 to form estimates of

{sk , lk , uk }4k=1 , and the results are shown in Table 4.

Table 4: Estimated power sums of eigenvalues for the DA chain
for Bayesian linear regression with non-Gaussian errors
k

Est. sk

√
Est. D/ N

Est. lk

Est. uk

1

35.587

0.121

0.000

34.587

2

2.465

0.020

0.042

1.210

3

1.325

0.014

0.222

0.687

4

1.102

0.012

0.313

0.564

The asymptotic 95% CIs for l4 and u4 are (0.241, 0.383) and (0.532, 0.597), respectively. Also, using
Bonferroni, we may conclude that asymptotically, with at least 95% confidence, λ1 ∈ (0.241, 0.597). The
24

Monte Carlo sample size required to secure a reasonable estimate does increase with n and p. For example,
in another example we considered where n = 20 and p = 4, we needed a Monte Carlo sample size of 1×107
to get decent results.

Acknowledgment. The second and third authors were supported by NSF Grant DMS-15-11945.

Appendix
A Proof of Theorem 1
Theorem 1. The DA operator P is trace-class if and only if
Z

SU

p(u, u) µ(du) < ∞.

(4)

If (4) holds, then for any positive integer k,
sk :=

κ
X

λki

=

Z

SU

i=0

p(k) (u, u) µ(du) < ∞.

(5)

2
Proof. Note that P is self-adjoint and non-negative. Let {gi }∞
i=0 be an orthonormal basis of L (πU ). The

operator P is defined to be trace-class if (see e.g. Conway, 2000)
∞
X
hP gi , gi iπU < ∞.

(25)

i=0

This condition is equivalent to P being compact with summable eigenvalues. To show that P being traceclass is equivalent to (4), we will prove a stronger result, namely
∞
X
i=0

hP gi , gi iπU =

Z

p(u, u)µ(du).

(26)

SU

We begin by defining two new Hilbert spaces. Let L2 (πV ) be the Hilbert space consisting of functions
that are square integrable with respect to the weight function πV (·). For f, g ∈ L2 (πV ), their inner product
is defined, as usual, by
hf, giπV =

Z

f (v)g(v)πV (v) ν(dv).

SV

25

Let L2 (πU × πV ) be the Hilbert space of functions on SU × SV that are square integrable with respect to the
weight function πU (·)πV (·). For f, g ∈ L2 (πU × πV ), their inner product is
hf, giπU ×πV =

Z

f (u, v)g(u, v)πU (u)πV (v) µ(du) ν(dv).

SU ×SV

2
Note that L2 (πV ) is separable. Let {hj }∞
j=0 be an orthonormal basis of L (πV ). It can be shown that

{gi hj }(i,j)∈Z2+ is an orthonormal basis of L2 (πU × πV ). Of course, gi hj denotes the function given by
(gi hj )(u, v) = gi (u)hj (v).
The inequality (4) is equivalent to


Z
πU,V (u, v) 2
πU (u)πV (v) µ(du) ν(dv) < ∞,
πU (u)πV (v)
SU ×SV
which holds if and only if the function ϕ : SU × SV → R given by
ϕ(u, v) =

πU,V (u, v)
πU (u)πV (v)

is in L2 (πU × πV ). Suppose (4) holds. Then by Parseval’s identity,
Z
p(u, u) µ(du) = hϕ, ϕiπU ×πV
SU

=

X

(i,j)∈Z2+

=

X

(i,j)∈Z2+

=

|hϕ, gi hj iπU ×πV |2
Z

gi (u)hj (v)πU,V (u, v) µ(du) ν(dv)

SU ×SV

∞ Z
∞ X
X
i=0 j=0

2

SV

Z

SU


2
gi (u)πU |V (u|v) µ(du) hj (v)πV (v) ν(dv) .

Again by Parseval’s identity, this time applied to the function on SV (and in fact, in L2 (πV ) by Jensen’s
inequality) given by
ϕi (v) =

Z

SU

gi (u)πU |V (u|v) µ(du),

26

we have
Z
∞
∞ X
X
|hϕi , hj iπV |2
p(u, u) µ(du) =
SU

i=0 j=0

=
=
=
=

∞
X

hϕi , ϕi iπV

i=0
∞ Z
X

i=0 SV
∞ Z
X
i=0 SV
∞ Z
X
i=0

SU

Z

2

SU

Z

SU

Z

gi (u)πU |V (u|v) µ(du) πV (v) ν(dv)

Z

SU

SU


gi (u′ )πU |V (u′ |v)πV |U (v|u) µ(du′ ) gi (u)πU (u) µ(du) ν(dv)


p(u, u′ )gi (u′ ) µ(du′ ) gi (u)πU (u) µ(du).

(27)

Note that the use of Fubini’s theorem in the last equality can be easily justified by noting that gi ∈ L2 (πU ),
P
and making use of Jensen’s inequality. But the right hand side of (27) is precisely ∞
i=0 hP gi , gi iπU . Hence,
R
(26) holds when SU p(u, u) µ(du) is finite.
To finish our proof of (26), we’ll show (25) implies (4). Assume that (25) holds. Tracing backwards

along (27) yields
X

(i,j)∈Z2+

|hϕi , hj iπV |2 < ∞.

This implies that the function
ϕ̃ :=

X

(i,j)∈Z2+

hϕi , hj iπV gi hj

is in L2 (πU × πV ). Recall that (4) is equivalent to ϕ being in L2 (πU × πV ). Hence, it suffices to show that
ϕ̃(u, v) = ϕ(u, v) almost everywhere. Define a linear transformation T : L2 (πU ) → L2 (πV ) by
Z
f (u)πU |V (u|v) µ(du), ∀f ∈ L2 (πU ).
T f (v) =
SU

By Jensen’s inequality, T is bounded, and thus, continuous. For any g =

27

P∞

i=0 αi gi

∈ L2 (πU ) and h =

P∞

j=0 βj hj

∈ L2 (πV ),
Z

SV

Z

ϕ(u, v)g(u)h(v) πU (u)πV (v) µ(du) ν(dv)
SU

= hT g, hiπV
=
=

∞
∞ X
X

i=0 j=0
∞ X
∞
X
i=0 j=0

αi βj hT gi , hj iπV
αi βj hϕi , hj iπV

= hϕ̃, ghiπU ×πV
Z Z
ϕ̃(u, v)g(u)h(v) πU (u)πV (v) µ(du) ν(dv),
=
SV

SU

where g ∈ L2 (πV ) is given by g(u) := g(u), and gi is defined similarly for i ∈ Z+ . This implies that for
any C1 ∈ U and C2 ∈ V,
Z
Z
ϕ(u, v) πU (u)πV (v) µ(du) ν(dv) =

ϕ̃(u, v) πU (u)πV (v) µ(du) ν(dv).

C1 ×C2

C1 ×C2

Note that
Z

1/2

SU ×SV

|ϕ̃(u, v)| πU (u)πV (v) µ(du) ν(dv) ≤ hϕ̃, ϕ̃iπU ×πV < ∞.

(28)

By (28) and the dominated convergence theorem, one can show that
Z
Z
o
n
ϕ̃(u, v) πU (u)πV (v) µ(du) ν(dv)
ϕ(u, v) πU (u)πV (v) µ(du) ν(dv) =
A := C ∈ U × V
C

C

is a λ system. An application of Dynkin’s π-λ theorem reveals that U × V ⊂ A. Therefore, ϕ̃(u, v) =
ϕ(u, v) almost everywhere, and (4) follows.
For the rest of the proof, assume that P is trace-class. This implies that P is compact, and thus admits
the spectral decomposition (see e.g. Helmberg, 2014, §28 Corollary 2.1) given by
Pf =

κ
X
i=0

λi hf, fi iπU fi ,

f ∈ L2 (πU )

where fi , i = 0, 1, . . . , κ, is the normalized eigenfunction corresponding to λi . By Parseval’s identity,
κ
∞ X
∞
X
X
λj |hgi , fj iπU |2
hP gi , gi iπU =
i=0 j=0

i=0

=
=

κ
X

j=0
κ
X
j=0

28

λj hfj , fj iπU
λj .

(29)

This equality is in fact a trivial case of Lidskii’s theorem (see e.g. Erdös, 1974; Gohberg et al., 2012). It
follows that (5) holds for k = 1.
We now consider the case where k ≥ 2. By (29) and a simple induction, we have the following decom-

position for P k .

k

P f=

κ
X
i=0

Hence P k

λki hf, fi iπU fi ,

f ∈ L2 (πU ) .

is trace-class with ordered positive eigenvalues {λki }κi=0 . Note that P k is a Markov operator whose

Mtd is p(k) (u, ·), u ∈ SU . Thus, in order to show that (5) holds for k ≥ 2, it suffices to verify P k is a DA

operator, for then we can treat P k as P and repeat our argument for the k = 1 case. To be specific, we’ll
show that there exists a random variable Ṽ taking values on SṼ , where (SṼ , Ṽ, ν̃) is a σ-finite measure space
and Ṽ is countably generated, such that for u ∈ SU ,
p

(k)

(u, ·) =

Z

SṼ

πU |Ṽ (·|v)πṼ |U (v|u) ν̃(dv),

(30)

where πṼ (·), πU |Ṽ (·|·), and πṼ |U (·|·) have the apparent meanings.
Let (Uk , Vk )∞
k=0 be a Markov chain. Suppose that U0 has pdf πU (·), and for any non-negative integer k,
let Vk |Uk = u have pdf πV |U (·|u), and let Uk+1 |Vk = v have pdf πU |V (·|v). It’s easy to see {Uk }∞
k=0 is a
stationary DA chain with Mtd p(u, ·). Suppose k is even. The pdf of Uk |U0 = u is
p

(k)

(u, ·) =

Z

SU

p(k/2) (u, u′ )p(k/2) (u′ , ·)µ(du).

Meanwhile, since the chain is reversible and starts from the stationary distribution, U0 |Uk/2 = u has the
same distribution as Uk/2 |U0 = u, which is just p(k/2) (u, ·). Thus, (30) holds with Ṽ = Uk/2 . A similar

argument shows that when k is odd, (30) holds with Ṽ = V(k−1)/2 .

B Proof of Proposition 6
Proposition 6. Suppose that h(·) is strictly positive in a neighborhood of the origin. If ω(z) can be written
Q
as ni=1 ωi (zi ), and there exists ξ ∈ (1, 4/3) such that for all i ∈ {1, 2, . . . , n},
Z

R+

u3d/2 h3 (u)
du < ∞,
R ξu
( 0 v d/2 h(v) dv)3 ωi2 (u)

then (15) holds, and thus by Theorem 2, second moment exists for the estimator (13).

29

Proof. Let Sd be the set of d × d positive definite matrices. For any β ∈ Rp , Σ ∈ Sd , z ∈ Rn , and
ξ ∈ (1, 4/3),
πU |V (β, Σ|z)πV3 |U (z|β, Σ)
Q
|Σ|−(n+d+1)/2 ni=1 exp{−zi (yi − β T xi )T Σ−1 (yi − β T xi )/2}
×
=R R
Qn
−(n+d+1)/2
T
T −1
T
i=1 exp{−zi (yi − β̃ xi ) Σ̃ (yi − β̃ xi )/2} dΣ̃ dβ̃
Rp Sd |Σ̃|
n
Y

3d/2

n
Y

zi

z
exp{−3zi (yi − β T xi )T Σ−1 (yi − β T xi )/2}h3 (zi )
R ∞i
{ 0 v d/2 exp[−v(yi − β T xi )T Σ−1 (yi − β T xi )/2]h(v) dv}3
i=1
Q
|Σ|−(n+d+1)/2 ni=1 exp{−zi (yi − β T xi )T [Σ/(4 − 3ξ)]−1 (yi − β T xi )/2}
≤R R
×
Qn
−(n+d+1)/2
T
T −1
T
i=1 exp{−zi (yi − β̃ xi ) Σ̃ (yi − β̃ xi )/2} dΣ̃ dβ̃
Rp Sd |Σ̃|
3d/2 3
h (zi )

.
R ξzi
d/2 h(v) dv)3
(
v
i=1
0

Note that
Z

−(n+d+1)/2

Sd

|Σ|

n
Y
i=1

−nd/2

= (4 − 3ξ)

Z

Sd

exp



T
zi
−
y i − β T xi
2

|Σ|−(n+d+1)/2

n
Y

exp

i=1

n



−

Σ
4 − 3ξ

−1

Rp

Z

Sd

πU |V (β, Σ|z)πV3 |U (z|β, Σ) dΣ dβ

−nd/2

≤ (4 − 3ξ)

The result follows immediately.

y i − β xi





dΣ

o
T
zi
yi − β T xi Σ−1 yi − β T xi dΣ.
2

Thus,
Z

T

n
Y
i=1

3d/2 3
h (zi )

zi

.
R ξz
( 0 i v d/2 h(v) dv)3

References
A HUES , M., L ARGILLIER , A. and L IMAYE , B. (2001). Spectral Computations for Bounded Operators.
CRC Press.
A LBERT, J. H. and C HIB , S. (1993). Bayesian analysis of binary and polychotomous response data. Journal
of the American statistical Association 88 669–679.
B ILLINGSLEY, P. (1995). Probability and Measure. 3rd ed. John Wiley & Sons.
B RISLAWN , C. (1988). Kernels of trace class operators. Proceedings of the American Mathematical Society
104 1181–1190.
C HAKRABORTY, S. and K HARE , K. (2017). Convergence properties of Gibbs samplers for Bayesian probit
regression with proper priors. Electronic Journal of Statistics 11 177–210.
30

C HAN , K. S. and G EYER , C. J. (1994). Discussion: Markov chains for exploring posterior distributions.
The Annals of Statistics 22 1747–1758.
C HOI , H. M. and H OBERT, J. P. (2013). Analysis of MCMC algorithms for Bayesian linear regression with
Laplace errors. Journal of Multivariate Analysis 117 32–40.
C HOI , H. M. and ROM ÁN , J. C. (2017). Analysis of Polya-Gamma Gibbs sampler for Bayesian logistic
analysis of variance. Electronic Journal of Statistics 11 326–337.
C ONWAY, J. B. (1990). A Course in Functional Analysis. 2nd ed. Springer-Verlag.
C ONWAY, J. B. (2000). A Course in Operator Theory. American Mathematical Soc.
D IACONIS , P., K HARE , K. and S ALOFF -C OSTE , L. (2008). Gibbs sampling, exponential families and
orthogonal polynomials (with discussion). Statistical Science 23 151–200.
D IACONIS , P. and S TROOCK , D. (1991). Geometric bounds for eigenvalues of Markov chains. The Annals
of Applied Probability 1 36–61.
E RD ÖS , J. (1974). On the trace of a trace class operator. Bulletin of the London Mathematical Society 6
47–50.
F ERNANDEZ , C. and S TEEL , M. F. (1999). Multivariate student-t regression models: Pitfalls and inference.
Biometrika 86 153–167.
G ARREN , S. T. and S MITH , R. L. (2000). Estimating the second largest eigenvalue of a Markov transition
matrix. Bernoulli 6 215–242.
G OHBERG , I., G OLDBERG , S. and K RUPNIK , N. (2012). Traces and Determinants of Linear Operators,
vol. 116. Birkhäuser.
H ELMBERG , G. (2014). Introduction to Spectral Theory in Hilbert Space. Elsevier.
H OBERT, J. P. (2011). The data augmentation algorithm: Theory and methodology. In Handbook of Markov
Chain Monte Carlo (S. Brooks, A. Gelman, G. Jones and X.-L. Meng, eds.). Chapman & Hall/CRC Press.
H OBERT, J. P. and M ARCHEV, D. (2008). A theoretical comparison of the data augmentation, marginal
augmentation and PX-DA algorithms. The Annals of Statistics 36 532–554.
K HARE , K. and H OBERT, J. P. (2011). A spectral analytic comparison of trace-class data augmentation
algorithms and their sandwich variants. The Annals of Statistics 39 2585–2606.
31

KOLTCHINSKII , V. and G IN É , E. (2000). Random matrix approximation of spectra of integral operators.
Bernoulli 6 113–167.
KONTOYIANNIS , I. and M EYN , S. P. (2012). Geometric ergodicity and the spectral gap of non-reversible
Markov chains. Probability Theory and Related Fields 154 327–339.
L AWLER , G. F. and S OKAL , A. D. (1988). Bounds on the l2 spectrum for Markov chains and Markov
processes: A generalization of Cheegers inequality. Transactions of the American Mathematical Society
309 557–580.
L IU , C. (1996). Bayesian robust multivariate linear regression with incomplete data. Journal of the American
Statistical Association 91 1219–1227.
L IU , J. S., W ONG , W. H. and KONG , A. (1994). Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and augmentation schemes. Biometrika 81 27–40.
L IU , J. S. and W U , Y. N. (1999). Parameter expansion for data augmentation. Journal of the American
Statistical Association 94 1264–1274.
M IRA , A. and G EYER , C. J. (1999). Ordering Monte Carlo Markov chains. School of Statistics, University
of Minnesota. technical report .
PAL , S., K HARE , K. and H OBERT, J. P. (2017). Trace class markov chains for bayesian inference with
generalized double pareto shrinkage priors. Scandinavian Journal of Statistics 44 307–323.
Q IN , Q. and H OBERT, J. P. (2016). Trace-class Monte Carlo Markov chains for Bayesian multivariate linear
regression with non-Gaussian errors. arXiv:1602.00136 .
ROBERTS , G. O. and ROSENTHAL , J. S. (1997). Geometric ergodicity and hybrid Markov chains. Electronic Communications in Probability 2 13–25.
ROY, V. and H OBERT, J. P. (2007). Convergence rates and asymptotic standard errors for markov chain
monte carlo algorithms for bayesian probit regression. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 69 607–623.
S INCLAIR , A. and J ERRUM , M. (1989). Approximate counting, uniform generation and rapidly mixing
Markov chains. Information and Computation 82 93–133.
TANNER , M. A. and W ONG , W. H. (1987). The calculation of posterior distributions by data augmentation
(with discussion). Journal of the American statistical Association 82 528–540.
32

VAN

DYK , D. A. and M ENG , X.-L. (2001). The art of data augmentation (with discussion). Journal of

Computational and Graphical Statistics 10 1–50.

33

