Fast Sampling for Strongly Rayleigh Measures with
Application to Determinantal Point Processes
Chengtao Li
Stefanie Jegelka
Suvrit Sra

ctli@mit.edu
stefje@csail.mit.edu
suvrit@mit.edu

Massachusetts Institute of Technology, Cambridge, MA 02139

arXiv:1607.03559v1 [cs.LG] 13 Jul 2016

Abstract
In this note we consider sampling from (non-homogeneous) strongly Rayleigh probability measures. As an important corollary, we obtain a fast mixing Markov Chain sampler for Determinantal
Point Processes.

1

Introduction

Probability distributions over combinatorial families of subsets are important to a variety of problems
in machine learning and related areas. Notable examples include discrete probabilistic models [11,
22, 26, 37, 40] for use in computer vision, computational biology, and Natural Language Processing;
combinatorial bandit learning [14]; model compression [34]; and low-rank matrix approximations [27].
Consequently, significant recent attention has been paid to sampling rapidly from certain structured
discrete distributions [21, 36], as well as from determinantal point processes [3, 27, 28] and progress on
sampling by optimization [16, 32].
Amongst these distributions, a widely used class is that of log-submodular measures. Formally, for
sets S, T ⊆ V, a log-submodular measure π : 2V → R+ satisfies the inequality
log π (S) + log π ( T ) ≥ log π (S ∪ T ) + log π (S ∩ T ).

(1.1)

Log-submodular measures are useful to several applications in machine learning and computer vision
[5, 26]; more generally, submodular functions are widely important across machine learning [4, 6, 25].
In this note, we focus on a specific subclass of log-submodular measures, namely, strongly Rayleigh
(SR) measures. These measures are intimately related to stable polynomials, a viewpoint first established
in [7], which has proved key to uncovering their remarkable properties, both for modeling as well as
for fast sampling. For instance, these measures exhibit negative association, a strong, “robust” notion of
negative dependence (we formally define SR measures in Section 2).
We mention below some important examples of SR measures.
Determinantal Point Processes. A Determinantal Point Process (DPP) is a measure over subsets given
by the principal minors of a positive semidefinite matrix K ∈ R N × N with eigenvalues in [0, 1]. Its
marginal probabilities satisfy
Pr(S ⊆ T ) = det(KS ),
(1.2)
where KS is the submatrix indexed by the elements in S, and T is the random set distributed as a DPP.
DPPs arise in random matrix theory, combinatorics, machine learning, matrix approximations, and
many other areas; see e.g., [8–10, 13, 23, 26, 28–31, 38].
(Weighted) regular and balanced matroids. The uniform distribution over the bases of certain matroids (regular matroids and balanced matroids [17, 35]) is SR, most notably, the uniform distribution
over spanning trees in a graph. Here, spanning trees are viewed as subsets of edges, and the distribution is over subsets of edges.

1

Product measures / Bernoullis conditioned on their sum. Assume there is a weight qi ∈ [0, 1] for each
element i ∈ V. The product measure π (S) = ∏i∈S qi ∏ j/
∈S (1 − q j ) is SR, as is its conditioning on sets of
a specific cardinality k, i.e., π 0 (S) = π (S | |S| = k) or π 0 (S) = 0 if |S| 6= k, and π 0 (S) ∝ π (S) otherwise.
Strongly Rayleigh measures have been underlying recent progress in approximation algorithms [1,
15, 20, 27], graph sparsification [18, 39], extensions to the Kadison-Singer problem [2], finite extensions
to free probability [33], and concentration of measure results [35].
Contributions. Despite their importance, efficient sampling methods are only known for special cases
of SR measures. In this note, we derive a provably fast mixing Markov Chain for efficiently sampling
general SR measures. For our analysis, we use the recent result of [3] (that analyzes fast mixing for
the subclass of k-homogeneous SR measures), along with the closedness properties of SR measures
established in the landmark work [7].

2 Sampling from Strongly Rayleigh Distributions
Strongly Rayleigh (SR) distributions capture the strongest form of negative dependence, while enjoying
a host of other notable properties [7]. Several important distributions exhibit the strong Rayleigh
property, for example, uniform distributions over spanning trees in graphs, and more generally, the
widely occurring Determinantal Point Processes. A distribution is strongly Rayleigh if its generating
polynomial pπ : C N → C,
p π ( z ) = ∑ π ( S ) ∏ zi
(2.1)
S ⊆V

i ∈S

is real stable. This means that if =(zi ) > 0 for all arguments zi of pπ (z), then pπ (z) 6= 0.
Markov Chain Sampling. We sample from π via a Markov Chain Monte Carlo method (MCMC),
i.e., we run a Markov Chain with state space 2V (the power set of V). All the chains discussed here are
ergodic. The mixing time of the chain indicates the number of iterations t that we must perform (after
starting from an arbitrary set S0 ∈ 2V ) before we can consider St a valid sample from π. Formally, if
δS0 (t) is the total variation distance between the distribution of St and π after t steps, then τS0 (ε) =
min{t : δS0 (t0 ) ≤ ε, ∀t0 ≥ t} is the mixing time to sample from a distribution e-close to π in terms of
total variation distance. We say that the chain mixes fast if τS0 is polynomial in N.
Existing samplers. Efficient sampling techniques have been studied for special cases of SR distributions. A popular method for sampling from Determinantal Point Processes uses the spectrum of the
defining kernel [23]. Generic MCMC samplers can also be derived, for example, previous work used
a simple add-delete Metropolis-Hasting chain [24]. Starting with an arbitrary set S ⊆ V, we sample a
point t ∈ V uniformly at random. If t ∈ S, we remove t with probability min{1, π (S \ {t})/π (S)}; if
t∈
/ S, we add it to S with probability min{1, π (S ∪ {t})/π (S)}. Algorithm 1 shows the (lazy) Markov
chain.
The add-delete chain can work well in practice [24], however, it does not always mix fast. An
elementary Determinantal Point Process has non-zero measure only on sets of a fixed cardinality; for
such a process (or a process close to it), the chain will stall or mix slowly.
Another special case of SR distributions are homogeneous SR measures. These measures are nonzero
only for some sets of a fixed cardinality k. Examples include Bernoulli distributions conditioned on
cardinality, uniform distributions on the bases of balanced matroids [17], and k-Determinantal Point
Processes. A natural MCMC sampler for these processes takes swapping steps: given a current set
S ⊆ V, it picks, uniformly at random, points s ∈ S and t ∈
/ S, and swaps them with probability
min{1, π (S ∪ {t} \ {s})/π (S)}. Algorithm 2 formalizes this procedure. Building upon results in [17],

2

Algorithm 1 Add/delete (Metropolis-Hasting) sampler
Require: SR distribution π
Initialize S ⊆ V
while not mixed do
Let b = 1 with probability 21
if b = 1 then
Pick t ∈ V uniformly at random
if t ∈ S then
S = S\{t} with probability min{1, π (S \ {t})/π (S)}
else
S = S ∪ {t} with probability min{1, π (S ∪ {t})/π (S)}
end if
else
Do nothing
end if
end while
Algorithm 2 Gibbs exchange sampler
Require: Homogeneous SR distribution π
Initialize S ⊆ V, π (S) > 0
while not mixed do
Let b = 1 with probability 21
if b = 1 then
Pick s ∈ S and t ∈
/ S uniformly randomly
S = S ∪ {t}\{s} with probability min{1, π (S ∪ {t} \ {s})/π (S)}
else
Do nothing
end if
end while

Anari et al. [3] recently showed that the mixing time for the swap sampler for homogeneous SR measures is polynomial in N, k, and log( eπ (1S ) ). These results are restricted to homogeneous SR measures,
0
and do not hold for arbitrary SR measures.

2.1

A fast mixing chain for general SR measures

In this note, we define a projection chain that works for arbitrary SR measures, and whose mixing time
is polynomial in N, k, and log( eπ (1S ) ). In particular, we make the results in [3, 17] accessible to general
0
SR measures by using specific closure properties [7].
The resulting Markov Chain is shown in Algorithm 3. Interestingly, this sampler uses a mixture of
add-delete and swap steps. Hence, intuitively, it preserves the good properties of either type of step.
In general the sampled sets can have arbitrary cardinality, and hence add-delete steps are needed. If
the distribution concentrates on a certain cardinality, the swap steps gain importance.
This intuition is supported by the following theorem.
Theorem 1. If π is a SR measure, the mixing time τS0 (ε) of the Markov chain in Algorithm 3 is given by




N
τS0 (ε) ≤ 2N 2 log
+ log(π (S0 ))−1 + log ε−1 .
(2.2)
| S0 |

3

Algorithm 3 Markov Chain for Strongly Rayleigh Distribution
Require: SR distribution π
Initialize R0 ⊆ [2N ] where | R0 | = N and take S = R0 ∩ V
while not mixed do
draw q ∼ Unif [0, 1]
draw t ∈ V \S and s ∈ S uniformly randomly
if q ∈ [0,

( N −|S|)2
)
2N 2

then

S = S ∪ {t} with probability min{1,
else if q ∈

−|S|)2 N −|S|
[ ( N2N
, 2N )
2

π (S∪{t})
π (S)

. Add t

π (S∪{t}\{s})
}
π (S)

. Exchange s with t

then

S = S ∪ {t}\{s} with probability min{1,
else if q ∈

|S|+1
}
N −|S|

×

−|S| |S|2 + N ( N −|S|)
,
)
[ N2N
2N 2

then

S = S\{s} with probability min{1,
else
Do nothing
end if
end while

π (S\{s})
π (S)

×

|S|
}
N −|S|+1

. Delete s

We may choose the initial set such that S0 makes the first term in the sum logarithmic in N (S0 =
R0 ∩ V in Algorithm 3).
Theorem 1 and Algorithm 3 make use of the closure of SR measures under symmetric homogenization
[7]. The idea underlying this construction is to introduce a “shadow” V 0 of the ground set V, and to
construct an N-homogeneous SR measure πsh on this joint ground set V ∪ V 0 . Importantly, the marginal
distribution on V under this joint measure is exactly π. The homogeneous measure πsh leads to a fast
mixing chain that is, however, not practical to implement. Hence, we reduce it to an equivalent, more
efficient chain.
Proof. First, we construct a symmetric homogenization of π, a measure πsh on V ∪ V 0 :
(
−1
if | R| = N;
π ( R ∩ [ N ])( R∩[NN ])
πsh ( R) =
0
otherwise.
If π is SR, so is its symmetric homogenization πsh . We use this property to derive a fast-mixing chain.
The results in [3] show that a Markov Chain with swap steps mixes rapidly for πsh . Precisely, they
show that for any k-homogeneous SR distribution on a ground set of size M, a Gibbs-exchange sampler
has mixing time
τR0 (ε) ≤ 2k( M − k)(log πsh ( R0 )−1 + log ε−1 ).
Here, M = 2N and k = N, leading to a mixing time of 2N 2 (log πsh ( T0 )−1 + log ε−1 ), or, equivalently,




N
τS0 (ε) ≤ 2N 2 log
+ log(π (S0 ))−1 + log ε−1 .
(2.3)
| S0 |
It remains to show that the chain in Algorithm 3 is equivalent to the Gibbs-exchange sampler for
πsh . In fact, one may be tempted to implement the exchange sampler directly. However, it doubles the
size of the ground set to 2N, and always maintains a set of size N. If N is large, this can be impractical.
Claim 1. The mixing time of Markov chain in Algorithm 3 has the same bound as Eq. (2.3).
1
2,

Our exchange sampler maintains a set R of cardinality | R| = N. In each iteration, with probability
the sampler does nothing, otherwise it proceeds. If it proceeds, it picks s ∈ R and t ∈ [2N ]\ R

4

Convergence for Inference

0.6
0.5

Potential Scale Reduction Factor

1.3

Marg
Cond-5
Cond-10

Potential Scale Reduction Factor

1.3

Add-Delete
Mix

1.25

1.25

1.2

1.2

Add-Delete
Mix

PSRF

Val

0.3
0.2
0.1

PSRF

0.4

1.15

1.15

1.1

1.1

1.05

1.05

0
-0.1

1

-0.2
0

0.5

1
# Iter

1.5
×10

1
0

2

2

4

6
# Iter

4

8

10
×10 4

(b)

(a)

0

2

4

6

8

# Iter

10
×10 4

(c)

Figure 1: (a) Convergence of marginal and conditional probabilities by Dpp on uniform matroid, (b,c)
comparison between add-delete chain (Algorithm 1) and projection chain (Algorithm 3) for two instances: slowly decaying spectrum and sharp step in the spectrum.
uniformly at random, and exchanges them with probability


π ( R ∪ {t} \ {s})
.
min 1, sh
πsh ( R)

(2.4)

If the exchange is accepted, then the new set is R ∪ {t} \ {s}.
To consider the projection of this chain onto V, let S = R ∩ V, and T = V \ R. There are in total four
possibilities for locations of s and t:
|S|( N −|S|)
, s ∈ S and t ∈ T, and we switch assignment of s and t with probability
N2
πsh ( R∪{t}\{s})
t}\{s})
} = min{1, π (S∪{
}. This is equivalent to switching elements between S
πsh ( R)
π (S)

1. With probability

min{1,
and T, i.e., an exchange step on V.
|S|( N −|S|)

2. With probability
, we have s ∈
/ S and t ∈
/ T. In this case, independent of whether we
N2
exchange s and t or not, the set S = R ∩ V remains the same. Hence, in this step, S remains
unchanged.
3. With probability
|S|
}.
N −|S|+1

we have s ∈ S and t ∈
/ T, and we switch with probability min{1,

π (S\{s})
π (S)

×

π (S∪{t})
π (S)

×

This is equivalent to deleting element s from S.

4. With probability
|S|+1
}.
N −|S|

| S |2
,
N2

( N −|S|)2
,
N2

we have s ∈
/ S and t ∈ T, and switch with probability min{1,

This is equivalent to adding element t to S.

Algorithm 3 performs those steps with exactly the same probabilities; hence, it is a projection of the
exchange chain for πsh and has the same mixing time.
Remarks. By using the SR property, we obtain a clean bound for fast mixing. In certain cases,
the above chain may mix slower in practice than a pure add-delete chain, since it is “lazier”, i.e., its
probability of stalling is higher. However, it is guaranteed to mix well, and, in other cases, can mix
much faster than the pure add-delete chain in [21, 24]. We observe both phenomena in our experiments.

3

Experiments

Next, we empirically study the how fast our samplers converge. We compare the strongly-Rayleigh
chain in Algorithm 3 (Mix) against a simple add-delete chain (Add-Delete). To monitor the convergence
5

of these Markov chains, we use potential scale reduction factor (PSRF) [12, 19] that runs several chains
in parallel and compares within-chain variances to between-chain variances. Typically, PSRF is greater
than 1 and will converge to 1 in the limit; if it is close to 1 we empirically conclude that chains have
mixed. Throughout experiments we run 10 chains in parallel for estimation, and declare “convergence”
at a PSRF of 1.05.
We use a Dpp on Ailerons data1 of size 200, and the corresponding PSRF is shown in Fig. 1b.
We observe that Mix converges slightly slower than Add-Delete since it is lazier. However, the AddDelete chain does not always mix fast. Fig. 1c illustrates a different setting, where we modify the
eigenspectrum of the kernel matrix: the first 100 eigenvalues are 500 and others 1/500. Such a kernel
corresponds to almost an elementary Dpp, where the size of the observed subsets sharply concentrates
around 100. Here, Add-Delete moves very slowly. Mix, in contrast, has the ability to exchange elements
and thus converges much faster than Add-Delete.

References
[1] Nima Anari and Shayan O. Gharan. Effective-resistance-reducing flows and asymmetric tsp. In
IEEE Symposium on Foundations of Computer Science (FOCS), 2015.
[2] Nima Anari and Shayan Oveis Gharan. The kadison-singer problem for strongly rayleigh measures and applications to asymmetric tsp. arXiv preprint arXiv:1412.1143, 2014.
[3] Nima Anari, Shayan Oveis Gharan, and Alireza Rezaei. Monte Carlo Markov chain algorithms
for sampling strongly Rayleigh distributions and determinantal point processes. COLT, 2016.
[4] Francis Bach. Learning with Submodular Functions: A Convex Optimization Perspective. Foundations
and Trends in Machine Learning, 2013.
[5] Olga Barinova, Victor Lempitsky, and Pushmeet Kohli. On detection of multiple object instances
using Hough transforms. IEEE Trans. on Pattern Analysis and Machine Intelligence, 9:1773–1784,
2012.
[6] Jeff Bilmes. Deep mathematical properties of submodularity with applications to machine learning. Tutorial at the Conference on Neural Information Processing Systems (NIPS), 2013.
[7] Julius Borcea, Petter Brändén, and Thomas Liggett. Negative dependence and the geometry of
polynomials. Journal of the American Mathematical Society, 22(2):521–567, 2009.
[8] Alexei Borodin. Determinantal point processes. arXiv:0911.1153, 2009.
[9] Alexei Borodin and Vadim Gorin. Lectures on integrable probability, 2012.
[10] Alexei Borodin and Grigori Olshanski. Distributions on partitions, point processes, and the hypergeometric kernel. Communications in Mathematical Physics, 211(2):335–358, 2000.
[11] Alexandre Bouchard-Côté and Michael I. Jordan. Variational inference over combinatorial spaces.
In NIPS, 2010.
[12] Stephen P Brooks and Andrew Gelman. General methods for monitoring convergence of iterative
simulations. Journal of computational and graphical statistics, pages 434–455, 1998.
[13] Alexander I Bufetov. Infinite determinantal measures. arXiv:1207.6793, 2012.
[14] Nicolo Cesa-Bianchi and Gabor Lugosi. Combinatorial bandits. In COLT, 2009.
[15] Amit Deshpande, Luis Rademacher, Santosh Vempala, and Grant Wang. Matrix approximation
and projective clustering via volume sampling. Theory of Computing, 2:225–247, 2006.
[16] Stefano Ermon, Carla P Gomes, Ashish Sabharwal, and Bart Selman. Embed and project: Discrete
sampling with universal hashing. In NIPS, pages 2085–2093, 2013.
[17] Tomás Feder and Milena Mihail. Balanced matroids. In STOC, pages 26–38, 1992.
[18] Alan Frieze, Navin Goyal, Luis Rademacher, and Santosh Vempala. Expanders via random spanning trees. SIAM Journal on Computing, 43(2):497–513, 2014.
1 http://www.dcc.fc.up.pt/657

~ltorgo/Regression/DataSets.html

6

[19] Andrew Gelman and Donald B Rubin. Inference from iterative simulation using multiple sequences. Statistical science, pages 457–472, 1992.
[20] Shayan O. Gharan, Amin Saberi, and Mohit Singh. A randomized rounding approach to the
Traveling Salesman Problem. In IEEE Symposium on Foundations of Computer Science (FOCS), 2011.
[21] Alkis Gotovos, Hamed Hassani, and Andreas Krause. Sampling from probabilistic submodular
models. In NIPS, pages 1936–1944, 2015.
[22] Dorothy M. Greig, Bruce T. Porteous, and Allan H. Seheult. Exact maximum a posteriori estimation for binary images. Journal of the Royal Statistical Society, 51(2), 1989.
[23] J Ben Hough, Manjunath Krishnapur, Yuval Peres, and Bálint Virág. Determinantal processes and
independence. Probability Surveys, 2006.
[24] Byungkon Kang. Fast determinantal point process sampling with application to clustering. In
NIPS, pages 2319–2327, 2013.
[25] Andreas Krause and Stefanie Jegelka. Submodularity in Machine Learning: New directions. Tutorial at the International Conference on Machine Learning (ICML), 2013.
[26] Alex Kulesza and Ben Taskar. Determinantal point processes for machine learning. Foundations and
Trends in Machine Learning. Now, 2012.
[27] Chengtao Li, Stefanie Jegelka, and Suvrit Sra. Fast DPP sampling for Nyström with application to
kernel methods. ICML, 2016.
[28] Chengtao Li, Suvrit Sra, and Stefanie Jegelka. Gaussian quadrature for matrix inverse forms with
applications. ICML, 2016.
[29] Russell Lyons. Determinantal probability measures. Publications Mathématiques de l’Institut des
Hautes Études Scientifiques, 98:167–212, 2003.
[30] Russell Lyons. Determinantal probability: basic properties and conjectures. arXiv:1406.2707, 2014.
[31] Odile Macchi. The coincidence approach to stochastic point processes. Advances in Applied Probability, pages 83–122, 1975.
[32] Chris J Maddison, Daniel Tarlow, and Tom Minka. A* sampling. In NIPS, 2014.
[33] Adam W Marcus. Polynomial convolutions and (finite) free probability, 2016.
[34] Zelda Mariet and Suvrit Sra. Diversity networks. In ICLR, 2016.
[35] Robin Pemantle and Yuval Peres. Concentration of lipschitz functionals of determinantal and
other strong rayleigh measures. Combinatorics, Probability and Computing, 23(01):140–160, 2014.
[36] Patrick Rebeschini and Amin Karbasi. Fast mixing for discrete point processes. COLT, 2015.
[37] David A. Smith and Jason Eisner. Dependency parsing by belief propagation. In EMNLP, 2008.
[38] Alexander Soshnikov. Determinantal random point fields. Russian Mathematical Surveys, 55(5):
923–975, 2000.
[39] Dan Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. In Symposium
on Theory of Computing (STOC), 2008.
[40] Jian Zhang, Josip Djolonga, and Andreas Krause. Higher-order inference for multi-class logsupermodular models. In ICCV, pages 1859–1867, 2015.

7

