Let µ : ( [n] k ) → R ≥0 be a density function on size k subsets of [n] = {1, . . . , n}, defining a distribution P[S] ∝ µ(S). The generating polynomial of µ is the multivariate k-homogeneous polynomial defined as follows:
We say that g µ is log-concave if log(g µ ) is a concave function over R n ≥0 . The study of log-concave polynomials has recently enabled breakthroughs on old conjectures about matroids, including the resolution of a conjecture of Mihail and Vazirani [MV89] on the expansion of the bases-exchange graphs [Ana+19], and Mason's ultra-log-concavity conjecture [Ana+18;BH18]. These results rely on the log-concavity of the generating polynomial for various distributions associated with matroids, most importantly the uniform distribution on the set of bases [AOV18].
Besides distributions associated with matroids, several other classes of distributions possess a log-concave generating polynomial. An important subclass consists of strongly Rayleigh distributions [BBL09] which includes determinantal point processes, distributions that have found numerous applications in machine learning [see KT12, for a survey]. A well-studied example belonging to all classes mentioned so far consists of the uniform distribution over spanning trees of a graph G = (V, E). Here n is the number of edges |E| in the graph and k is the number of edges in a spanning tree, i.e., |V| -1. Spanning trees of a graph form bases of a matroid called the graphic matroid [see, e.g., Oxl06] and they can also be viewed as a determinantal point process because of the matrix-tree theorem [see, e.g., BBL09], and are consequently strongly Rayleigh.
The motivation behind the conjecture of Mihail and Vazirani [MV89] was to solve the problem of approximately sampling from bases of a matroid. After this conjecture was made, efficient sampling algorithms were developed for various special classes of matroids [FM92; Gam99; JS02; Jer+04; Jer06; Clo10; CTY15; GJ18] until Anari, Liu, Oveis Gharan, and Vinzant [Ana+19] showed an efficient approximate sampling algorithm for all matroids. This algorithm used a variant of random walks on the so-called "basesexchange" graphs of matroids, that is known as the "down-up" random walk studied in the context of high-dimensional expanders [KM16; DK17; KO20]. For a distribution defined by µ : ( [n] k ) → R ≥0 , the down-up random walk P starts from a set S 0 ∈ ( [n] k ) and produces the Markovian sequence S 0 , S 1 , S 2 , . . . as follows:
Algorithm 1: Down-up Walk for t = 0, 1, 2, . . . do Let T t ∈ ( S t k-1 ) be a subset of S t obtained by dropping one element of S t uniformly at random. Let S t+1 = T t ∪ {e}, where the element e is chosen with probability ∝ µ(T t ∪ {e}).
The random walk P has µ as its stationary distribution and can be efficiently implemented by probing µ on at most n different sets each time. Thus assuming oracle access to µ, or in the case of matroids, an independence oracle for the matroid, each step of P takes O(n) time. The challenging part has been establishing the mixing time of P, i.e., bounds on the time t such that the distribution of S t is ϵ-close in total variation distance to the one defined by µ: t mix (P, S 0 , ϵ) := min {︁ t ⃓ ⃓ ∥P t (S 0 , •)µ(•)∥ TV ≤ ϵ }︁ .
Anari, Liu, Oveis Gharan, and Vinzant [Ana+19] proved that when µ has a log-concave generating polynomial, the spectral gap of the random walk P is at least 1/k. This implied that
Later, Cryan, Guo, and Mousa [CGM19] proved a Modified Log-Sobolev Inequality (MLSI) for the same random walk which resulted in a tighter mixing time:
These results lead to efficient algorithms assuming that the mass of the starting set, P µ [S 0 ], is not terribly small; this can often be achieved in practice. For example, for matroids, any starting basis S 0 will satisfy P µ [S 0 ] ≥ 1/( n k ) ≥ n -k , because the number of bases is at most ( n k ). Consequently the above bounds turn into t mix (P, S 0 , ϵ) ≤ O(k(k log(n) + log(1/ϵ))) and t mix (P, S 0 , ϵ) ≤ O(k(log k + log log n + log(1/ϵ))) respectively. However, for other distributions µ with a log-concave generating polynomial, even in the very special case of determinantal point processes, there is no control on min {P µ [S 0 ] | S 0 ∈ supp(µ)}, so one has to rely on clever tricks to find a good starting set S 0 ; even then, the best hope is to find a set S 0 with P µ [S 0 ] ≳ 1/( n k ), which results in a mixing time mildly depending on n. Historically, earlier works on a subclass of matroids, called balanced matroids, followed a similar development, where initially a spectral gap result was proved, resulting in a running time 1 of O(nk(k log n + log(1/ϵ))) followed by MLSI which resulted in a mixing time of O(k(log k + log log n + log(1/ϵ))) [see MT06, for a survey]. Noting that the term log log n seems unnecessary, Montenegro and Tetali [MT06] raised the question of proving a better inequality that would result in a running time of O(nk log(k/ϵ)). They specifically hoped for the possibility of proving a Nash inequality, an advanced type of functional inequality used to derive very tight mixing times for some Markov chains [MT06]. We believe there are barriers to using functional inequalities in general to prove O(k log(k/ϵ)) mixing time for the down-up random walk; we defer an explanation of this to a future version of this paper. However, without proving new functional inequalities, we manage to sidestep this barrier and improve the running time to the conjectured O(nk log(k/ϵ)) for not just balanced matroids, but the class of all matroids.
Our first result is a tight analysis of the mixing time, entirely removing the dependence on P µ [S 0 ] and n.
Theorem 1. For any distribution defined by µ : ( [n] k ) → R ≥0 with a log-concave generating polynomial g µ , the mixing time of the down-up random walk P, starting from any S 0 in the support of µ is
Note that generally we cannot hope for a better mixing time than k log k; each step of the random walk P replaces one element of the current set, and by a coupon collector argument, at least ≃ k log k steps are needed to replace every element of the starting set S 0 . As long as k is not too close to n, say k < 0.99n, replacing every starting element is needed for sufficient mixing, even for the simple distribution µ which is uniform over ( [n] k ). Our mixing time bound is an asymptotic improvement over prior work for k = O(1), or more generally when k is smaller than log(n) ϵ for all ϵ > 0. Another consequence of the new mixing time bound is that it enables the analysis of the down-up random walk when n is infinitely large; for example, this is the case for continuous determinantal point processes [see, e.g., OR18]. 2 To avoid complicating the notation, we do not consider infinitely large ground sets in this paper, but note that the results do generalize to such cases.
Our next result is the first quasi-linear time algorithm to sample from the uniform distribution over forests of a graph G = G(V, E). This improves upon the recent result by Anari, Liu, Oveis Gharan, and Vinzant [Ana+19] which gives a polynomial time algorithm to sample random forest, but the run-time of this algorithm is far from being linear in the number of edges. Their algorithm samples and counts forests of fixed-size k for each k ≤ |V| -1, thus takes at least Ω(|V| |E|) time. Moreover, they employ the approximate sampling to approximate counting reduction [JVV86; AD20], which introduces large polynomial blow-up in run-time. For applications of sampling random forests, see e.g. [Goe+].
In addition, we show a similar algorithm that also runs in quasi-linear time and samples from the uniform distribution over spanning trees of G. Much attention has been paid to the problem of sampling random spanning trees over the years, starting from the seminal works of Aldous [Ald90] and Broder [Bro89] who proposed a simple routine to extract a random spanning tree from the trace of a random walk on G itself. Subsequent works introduced improved algorithms [Wil96; CMN96; KM09; MST14; Dur+17a; Dur+17b] until finally Schild [Sch18] managed to obtain an almost-linear time algorithm running in time n 1+o(1) on graphs with n edges. This algorithm and that of several prior works were all based on the original work of Aldous [Ald90] and Broder [Bro89]; they achieved an improved running time by employing several clever, but complicated, tricks to shortcut the trace of a random walk over G. Our algorithms to sample a random spanning tree or random forest is wholly different, based on the down-up random walk, that achieves a nearly-linear running time of n log 2 (n), while being arguably much simpler to describe and implement. Our algorithms can be naturally extended to sample from weighted distributions over forests or spanning trees as well.
Theorem 2. There is an algorithm that takes a weighted graph G = G(V, E) on n edges with weight function w : E → R ≥0 , parameters q ≥ 0 and ϵ > 0 as input and outputs a forest F ⊆ E in time O(n log(n) log(n/ϵ)); the distribution of F is guaranteed to be ϵ-close in total variation distance to the distribution µ over forests of G defined by
where k is the rank of the graphic matroid of G, and |F| denotes the number of edges in F.
In particular, when w(e) = 1 for all e, µ is the uniform distribution on forests of G if q = 1, and is the uniform distribution on spanning trees of G if G is connected and q = 0.
In fact, we can extend Theorem 2 to allow sampling from the uniform distribution over forests of size between k 1 and k 2 , for any parameters k 1 , k 2 , in quasi-linear time.
Theorem 3. There is an algorithm that takes a weighted graph G = G(V, E) on n edges with weight function w : E → R ≥0 , parameters q ≥ 0, k 1 , k 2 ∈ N and ϵ > 0 as input and outputs a forest F ⊆ E in time O(n log(n) log(n/ϵ)); the distribution of F is guaranteed to be ϵ-close in total variation distance to the distribution
Since our algorithms are based on the MCMC method, they can only approximately sample from the forest or spanning tree distribution. In contrast, some of the prior works, including [Sch18], can sample exactly from the spanning tree distribution. This is mostly an inconsequential difference in practice, as no polynomial-time user of the algorithm can sense a difference between exact sampling and approximate sampling; one simply needs to set ϵ to be inverse-polynomially small.
We remark that our technique also leads to algorithm(s) that perform the more general task of approximately sampling from the uniform distribution over the family of independent sets of an arbitrary matroid, given access to suitable oracles. Specifically, for a matroid M = ([n], I) of rank k, an algorithm similar to the one from Theorem 2 samples from a distribution that is ϵ-close to the uniform distribution over the independent set I using O(n log n ϵ ) calls to a data structure O ′ that maintains a set S ⊆ [n], guaranteed to contain at most one circuit and allows for • Addition, and removal of an element from S, provided we maintain the property that S contains at most one circuit.
• Outputing a uniformly random element from the unique circuit in S if, such a circuit exists.
For graphic matroids, we can implement O ′ with amortized quasi-constant query time using link-cut trees. In general, since the input S is guaranteed to have size at most k + 1, we can implement each call to O ′ using O(k) calls to the more familiar independence oracle O I for M, resulting in a O(kn log n ϵ )-time algorithm.
In order to prove Theorem 1, our strategy is to combine a new analysis of the initial steps of the down-up random walk with the previously known Modified Log-Sobolev Inequality [CGM19]. Specifically we show that conditioned on having replaced every element of the starting set S 0 at least once by time t, the set at time t can be used as a warm start for the rest of the steps. Specifically, we show that the density of the set at time t w.r.t. µ, conditioned on this event, is upper-bounded by only a function of k.
In order to prove this, we introduce a new property of functions µ : ( [n] k ) → R ≥0 that we call α-approximate exchange. This property says that for every S, T ∈ ( [n] k ), and i ∈ S, there exists j ∈ T such that
Note that when µ takes values in {0, 1} and α ≥ 1, this property becomes equivalent to the famous strong basis exchange axiom of matroids [Oxl06]; if B = µ -1 (1) is the family of sets indicated by µ, this property says that for every S, T ∈ B and i ∈ S, there exists j ∈ T such that Si + j ∈ B and T + ij ∈ B. This property can be seen as a quantitative variant of strong basis exchange. Alternatively, it can be viewed as an approximate and multiplicative form of M ♮ -concavity, a cornerstone of discrete convex analysis [MS99].
We prove that every µ with a log-concave generating polynomial satisfies 2 O(k) -approximate exchange. Crucially, our α does not depend on n. We remark that Brändén and Huh [BH18] showed a result that can be thought of as a partial converse to this. They proved that M ♮ -concavity of log µ, equivalent to 1-approximate exchange property, implies that the generating polynomial of µ is log-concave.
We show that a similar approximate exchange property implies that a simple local search algorithm gives a k 2k -approximation on the problem of maximizing µ(S) for µ : ( [n] k ) → R ≥0 generated by a log-concave polynomial. If the generating polynomial of µ : ( [n] k ) → R ≥0 is moreover strongly Rayleigh, then µ satisfies a slightly stronger exchange property (see Lemma 24) that in turn implies greedy gives a k 2kapproximation of max S µ(S) (see Lemma 25). This is a generalization of Khachiyan [Kha95]'s classical result that greedy produces a k O(k) -approximation of the (sub)determinant maximization problem [Kha95; Di +14; Nik15], as well as [Pac04; ÇM10]'s more recent result that greedy gives k O(k) -approximation for the largest j-dimensional simplex problem. The best result on the largest j-dimensional simplex problem is a 2 O(k) -approximation by Nikolov [Nik15], matching the lower bound given by [Di +14; ÇM10].
We discuss the high-level ideas for proving Theorem 2. For the sake of simplicity, we consider the unweighted case, i.e., w(e) = 1 for all e ∈ E. It would be illuminating to first discuss the special case where q = 0, G is connected, and µ is uniformly distributed over spanning trees of G. We would like to use the down-up random walk from Theorem 1 to sample from µ. Though the down-up walk on the support of µ mixes in nearly-linear time, we do not see a way to implement each step of it in polylogarithmic time. Fortunately, the down-up random walk on an equivalent family of sets, the dual of the graphic matroid of G, which consists of the complements of spanning trees, also mixes fast, and we can implement each step in amortized O(log n)-time using the link-cut tree data structure [ST83; RTF18].
For q ̸ = 0, the distribution µ over forests of G is not homogeneous, i.e., the support of µ contains different-size subsets of E, so we cannot immediately apply Theorem 1. Let µ ¯be the complement distribution of µ i.e. µ ¯(E \ F) = µ(F) if F is a forest, then sampling from µ and from µ ¯are equivalent. We add auxiliary elements to each F ∈ supp(µ ¯) to obtain a homogeneous distribution. More precisely, we design a homogeneous distribution µ ↑ :
where Y is the set of auxiliary elements, such that the generating polynomial of µ ↑ is log-concave. Specifically, in Lemma 26, we prove, using results of Brändén and Huh [BH19], that for any matroid M of rank r over ground set
\S is strongly log-concave, and then use polarization (see Proposition 16) to transform f M into a multi-affine homogeneous log-concave polynomial
The fact that f M is strongly log-concave could be of independent interest.
The distribution µ ↑ is generated by f ↑ M . Our algorithm runs the down-up random walk on µ ↑ , which mixes rapidly by Theorem 1, and then outputs E \ (T t ∩ E) where T t ∈ supp(µ ↑ ) is the random set we obtained after t = O(n log n ϵ ) down-up steps. Each step of the walk, even in the weighted case, can again be implemented in amortized O(log n)-time using link-cut trees [ST83; RTF18].
If we only consider the effect of the down-up walk on T t,E := E \ (T t ∩ E), then each step of the down-up walk can be viewed as follows:
• With probability 1 -|T t,E | n , sample an edge e ̸ ∈ T t,E uniformly at random and add e to T t,E .
• If there is a cycle formed in T t,E by the previous operation, remove an edge uniformly at random from the cycle. Else, with probability q 1+q , remove an edge uniformly at random from T t,E . Note that this has no effects if T t,E is already empty.
Observe that if q = 0, we never remove an edge from T t,E unless T t,E contains a cycle, thus if T 0,E is a spanning tree then so is T t,E for all t. For q = 0, our algorithm (to sample random spanning trees) is the same as the one proposed by Russo, Teixeira, and Francisco [RTF18]. Despite not having the tight mixing time analysis, they empirically observed fast mixing times for the proposed algorithm, and additionally showed how link-cut trees can be used to implement each step.
To prove Theorem 3, we only need to show that the following polynomial is strongly log-concave:
Then, we apply this for the matroid M whose bases are the size-k 2 forests of graph G. Next, we employ the polarization trick then running down-up walk framework which we use to prove Theorem 2.
In Section 2 we provide some background on Markov chains and geometry of polynomials. In Section 3 we prove Theorem 1. In Section 4 we prove certain approximate exchange properties, and their algorithmic implications. In Section 5 we prove Theorems 2 and 3, as well as other results on sampling independent sets of an arbitrary matroid. The results in Section 5 can be read independently of Sections 3 and 4.
We use [n] to denote the set {1, . . . , n} and ( [n] k ) to denote the family of size k subsets of [n]. When n is clear from context, we use 1 S ∈ R n to denote the indicator vector of the set S ⊆ [n], having a coordinate of 0 everywhere except for elements of S, where the coordinate is 1. We use conv to denote the operator that maps a set of points to their convex hull.
We use z S as shorthand for {z i | i ∈ S} and z S as shorthand for ∏ i∈S z i . For polynomial f = ∑ S c S z S ∈ C[z 0 , • • • , z n ], we let the support of f be supp( f ) := {S : c S ̸ = 0}, and write ∂ i f as shorthand for ∂ f ∂z i .
We use e k (z 1 , • • • , z n ) to denote the k-th symmetric polynomial in z 1 , • • • , z n . We sometimes abuse notation and write e k (u, z S ) to denote the k-th symmetric polynomial in variables z S ∪ {u} .
In this paper we use one of the many cryptomorphic definitions of a matroid in terms of the polytope of its bases. For equivalence to other prominent definitions of a matroid, and more generally references to facts stated here see [Oxl06]. Definition 4. We say that a family B ⊆ ( [n] k ) is the family of bases of a matroid if the polytope conv {1 B | B ∈ B} has only edges of the minimum possible length, namely √ 2. We call k the rank of the matroid, and [n] the ground set of the matroid.
We let the family of independent sets of the matroid be I =
A well-known fact about matroids, that can be easily derived from Definition 4, is that the dual of a matroid, defined below, is another matroid.
is the family of bases of a matroid, then the following is also the family of bases of another matroid, called the dual matroid:
In this paper we will use a famous class of matroids constructed from graphs, called graphic matroids.
Proposition 6. Let G = (V, E) be a graph. Then the following is the family of bases of a matroid, called the graphic matroid of G:
{T ⊆ E | T forms a spanning forest}.
Note that the rank of the graphic matroid is ≤ |V| -1 and the ground set is E. If G is connected, then the bases are spanning trees of G, and the rank is exactly |V| -1.
Definition 7 (M-convex sets). We define a subset J ⊆ N n to be M-convex if it satisfies any one of the following equivalent conditions:
• For any α, β ∈ J and any index i satisfying α i > β i , there is an index j satisfying α j < β j and αe i + e j ∈ J.
• For any α, β ∈ J and any index i satisfying α i > β i , there is an index j satisfying α j < β j and αe i + e j ∈ J and βe j + e i ∈ J.
We note that any M-convex set J must be a subset of
Definition 8 (Half-plane stable). Consider an open half-plane
We say a polynomial g(z
In particular, the zero polynomial is H θ -stable.
We call H 0 and H π/2 the upper-half and right-half plane respectively. We say g is Hurwitz stable if it is H π/2 -stable. We say g is real stable if it is H 0 -stable and has real coefficients.
We observe that for homogeneous polynomials, the definition of H θ -stable is equivalent for all angles θ.
A distribution µ : 2 [n] → R ≥0 is strongly Rayleigh if and only if its generating polynomial is real stable [see BBL09]. Real stability is preserved under differentiation, and identification:
is real stable, then the following are also real stable
We will need the following classical fact [see e.g. BBL09, Proposition 3.1 for a proof].
Theorem 10. For k ≤ n, the k-th symmetric polynomial in n variables e k (z 1 , • • • , z n ) is real stable.
For a distribution or density function µ : ( [n] k ) → R ≥0 we denote by g µ the generating polynomial of µ defined as
We call a polynomial g ∈ R[z 1 , . . . , z n ] with nonnegative coefficients log-concave when viewed as a function, it is log-concave over the positive orthant, i.e., for x, y ∈ R n ≥0 and λ ∈ (0, 1)
For a multiaffine polynomial g, its derivatives can be obtained as
This shows that the derivatives of a multiaffine log-concave polynomial are limits of log-concave polynomials, which themselves are log-concave. It follows that a multiaffine homogeneous log-concave polynomial satisfies the seemingly stronger notions of strong log-concavity [Gur09] and complete log-concavity [AOV18; Ana+18; BH19]. The latter means that such polynomials are also closed under directional derivatives Lemma 11. Let g ∈ R[z 1 , . . . , z n ] be a multiaffine homogeneous polynomial with nonnegative coefficients. If g is log-concave, then it is completely log-concave as well, which means that for any k ∈ Z ≥0 and directions v 1 , . . . , v k ∈ R n ≥0 , the following polynomial is log-concave:
We will use an alternative characterization of strongly/completely log-concave homogeneous polynomial
In particular, if g is strongly log-concave then its support supp(g) is M-convex, and all its (directional) derivatives are strongly log-concave as well. For the equivalence of strongly log-concave, completely log-concave, and Lorentzian homoegeneous polynomials, see [BH19].
An important class of polynomials are those associated with uniform distributions over bases of a matroid.
k ) is the family of bases of a matroid, then the following polynomial is log-concave:
Theorem 14 ([BH19; Ana+18]). For any matroid M with family of independent sets I, the polynomial
We will use the following simple fact about log-concave polynomials:
Proposition 15 ([see Ana+18; BH18]). If g is a log-concave polynomial with nonnegative coefficients, then ∇ 2 g evaluated at any point in the positive orthant has at most one positive eigenvalue.
One of the basic operations preserving (complete) log-concavity is composition with a linear map. That is if T : R m → R n is an affine linear map for which T(R m ≥0 ) ⊆ R n ≥0 , then g • T is (completely) log-concave as well. We state other operations that preserve strong log-concavity/Lorentzian property.
The polarization map ∏ ↑ κ is a linear map that sends monomial z α = ∏ n i=1 z
where e α i ({z ij } 1≤j≤κ i ) is the elementary symmetric polynomial of degree α i in the variables {z ij } 1≤j≤κ i and
1≤i≤n is strongly log-concave then ∏ ↑ κ (g) is also strongly log-concave. Proposition 17. The product of two strongly log-concave polynomials is also strongly log-concave.
For two distributions ν, µ we define the Kullback-Leibler divergence, KL-divergence for short, D KL (ν ∥ µ) as
]︃ , and the total variation distance between ν and µ as ∥ν -
The two are related by Pinsker's inequality:
Proposition 18 ([see, e.g., CT12]). KL-divergence and the total variation distance are related by the following inequality
Cryan, Guo, and Mousa [CGM19] proved shrinkage of the KL-divergence under the down-up random walk. Coupled with Pinsker's inequality, this resulted in a mixing time bound.
k ) → R ≥0 are distributions where µ has a log-concave generating polynomial, and P is the down-up random walk operator whose stationary distribution is µ, then
In this section we prove Theorem 1 by analyzing the down-up random walk for distributions µ that have a log-concave generating polynomial. As a reminder, in each step, the down-up random walk transitions from a set S ∈ ( [n] k ) to S ′ ∈ ( [n] k ) as follows: • From S choose a subset T ⊆ S of size k -1 uniformly at random.
• From all supersets S ′ ⊇ T, choose one with probability ∝ µ(S ′ ).
Notice that the first step above simply drops a uniformly random element, and the second step replaces it with a new one (potentially the same element). Our high-level strategy is to prove that in O(k log k) steps, every element of the initial set is replaced at least once, and when this happens the distribution becomes a warm start and converges to µ in an additional O(k log k) steps.
Let τ be the first time such that every element in our initial set has been replaced at least once. In other words think of initial elements as unmarked, and every time we replace an element we mark the new element brought in. Then τ is the first time that every element is marked.
We will prove the following: Lemma 20. Let S t be the set at time t in the down-up random walk. Then for any X ∈ ( [n] k ) and any time t,
Note that without 2 O(k 2 ) , the r.h.s. is simply the stationary distribution. So this statement can be understood to say that as long as we have replaced each element at least once, we cannot be too far off from the stationary distribution.
Before proving Lemma 20, let us finish the proof of Theorem 1 assuming it.
Proof of Theorem 1 assuming Lemma 20. Note that for any fixed time t, we can simply bound P[τ > t] by k(1 -1/k) t ≤ ke -t/k . In particular this probability rapidly converges to 0 after about k log k steps. Now let t 1 < t 2 be two time indices. Let ν t denote the distribution of the state of random walk, i.e., S t , at time t. Our goal is to bound ∥ν t -µ∥ TV , where for simplicity of notation, we assume µ is properly normalized to be a probability distribution. Let ν ′ t be the distribution of S t conditioned on τ ≤ t, and let ν ′′ t be the distribution of S t conditioned on τ ≥ t. Then we can write
If P denotes the random walk operator, then note that ν t 2 = ν t 1 P t 2 -t 1 . So we get
Using the triangle inequality we can bound
Here we used the fact that P[τ ≤ t 1 ] ≤ 1, and ∥ν ′′ t 1 P t 2 -t 1 -µ∥ TV ≤ 1; the latter inequality is because ∥•∥ TV is always upper bounded by 1.
We can bound the second term in the above inequality by ke -t 1 /k as stated before. For the first term, note that the KL-divergence between ν ′ t 1 and µ is at most O(k 2 ) by Lemma 20. This is because
So by Lemma 19 in t 2t 1 steps this KL-divergence decreases to (1
. By Pinsker's inequality, Proposition 18, we get that
So in the end we get the following bound
In order for this to be at most ϵ, it is enough to make sure that min {t 1 , t 2t 1 } = Ω(k log k + k log 1 ϵ ). So we can simply let t 1 = t 2 /2, and then make sure that t 2 = Ω(k log(k/ϵ)).
As the main tool we use to prove Lemma 20, we introduce a new inequality for log-concave polynomials, that we call approximate exchange. We state the inequality below and defer its proof to Section 4.
k ) → R ≥0 with a log-concave generating polynomial satisfies a 2 O(k) -exchange property. That is, for every S, T ∈ ( [n] k ) and i ∈ S there exists j ∈ T such that
Armed with Lemma 21, let us prove Lemma 20.
Proof of Lemma 20. Let's look at the down-up walk process with orders. This means that we start with some elements e 1 , . . . , e k that together form the starting set. In each time step we replace one of the e i 's. But we keep track of the ordering and do not convert these to sets. So we can talk about e t i as the i-th element at time t. In particular S t is simply the unordered collection {e t 1 , . . . , e t k }. Let's say that X = { f 1 , . . . , f k }. Then to have S t = X, there must be some permutation of f 1 , . . . , f k that equals e t 1 , . . . , e t k . We will show that for any such permutation the promised bound in Lemma 20 holds. Since there are k! = 2 O(k log k) many permutations, this extra factor of k! can be absorbed into the factor of 2 O(k 2 ) without any loss. So we fix an arbitrary permutation, w.l.o.g. the identity permutation, and try to bound the following
Since we are conditioning on τ ≤ t, note that there must be some time τ i ≤ t, which is the last time before t where the i-th element gets replaced by the down-up random walk. We will bound the above probability, even conditioned on τ 1 , . . . , τ k having any set of fixed values up to t. Note that the index of the element that gets replaced in every step is uniformly random and independent of everything else that happens in the random walk, in particular the identity of the elements that come in as replacements. In the rest of the proof, we condition on the indices of the elements that get replaced at every step up to time t; note that this also uniquely determines τ 1 , . . . , τ k , so we assume τ 1 , . . . , τ k are some fixed time indices. W.l.o.g. assume that τ 1 < τ 2 < • • • < τ k . We will use induction to prove the following statement for i = 0, . . . , k:
Notice that for i = 0, both sides are trivially equal to 1, and for i = k, this inequality is the main statement we want to prove.
It remains to show the inductive step. We will show that going from i -1 to i, the l.h.s. gets multiplied by a smaller quantity compared to the r.h.s. If we have below inequality in hand, then it is not hard to see that we can complete the induction, since the factors that get multiplied on each side are the two sides of this inequality.
Instead of conditioning only on f 1 , . . . , f i-1 being chosen at the appropriate times on the l.h.s., we will refine the conditioning and condition on the history of the random walk up to time τ i -1. This means we can in particular assume that the elements e τ i i+1 , . . . , e τ i k are fixed, that e
, and the only uncertain thing is what the i-th element is being replaced by at time τ i .
Let S = { f 1 , . . . , f i , e τ i i+1 , . . . , e τ i k }. Then the conditional probability of choosing f i at time τ i is:
On the other hand
So we have to show the following:
We will give an injection from the terms on the l.h.s. to the terms in the expanded form of the r.h.s. Choose some set T ∋ f 1 , . . . , f i-1 . Apply Lemma 21 to S and T with the element f i ∈ S. We get that there must be some element e ∈ T such that
Note that V := Sf i + e contains Sf i , and
appears on the r.h.s. of the desired inequality. So for each T appearing on the l.h.s. of the desired inequality we produced a pair of U and V. Note that this mapping from T to (X, Y) is injective. This is because given (X, Y), we can recover T as the xor/symmetric difference of the other three sets, that is T = S∆U∆V.
In this section we prove Lemma 21. In addition, we show other variant(s) of approximate exchange property (Lemma 24) with implication for the approximation guarantee of local search and greedy on the problem of max S µ(S) given a log-concave/strongly-Rayleigh distribution µ :
Definition 22. We say that µ : ( [n] k ) → R ≥0 has an α-approximate exchange property, or α-exchange for short, if for every S, T ∈ ( [n] k ) and every i ∈ S, there exists j ∈ T such that
Note that if µ is the indicator of bases of a matroid, then it has a 1-exchange property, also known as the strong basis exchange property [Oxl06].
Although we do not directly need it, we give another example where approximate exchange can be proven by elementary means. This is the class of k-determinantal point processes [BBL09; KT12].
Proposition 23. Suppose that µ :
for some vectors v 1 , . . . , v n ∈ R k . Then µ has a k 2 -exchange property.
Proof. It is enough to consider the case where S and T are disjoint; otherwise, the problem can be reduced to lower values of k by taking out the intersection, and projecting all vectors on the orthogonal complement of the space spanned by the intersection.
Define the number β j as √︁ µ(Si + j)µ(T + ij) and let α be √︁ µ(S)µ(T). The Plücker relations for the Grassmanian [see, e.g., Abe80] say that a signed sum of α and β j is zero:
This means that there is at least one j such that |β j | ≥ 1 k α, and this concludes the proof.
Next we take steps to prove Lemma 21, namely that if µ : ( [n] k ) → R ≥0 has a log-concave generating polynomial g µ , then µ has a 2 O(k) -exchange property. We conjecture that a k O(1) -exchange property should hold, but even if true, this will not improve the mixing time results in this paper beyond constants hidden in the O(•) notation.
Our strategy is to prove the case of k = 2 of Lemma 21 by using log-concavity of g µ (note that k = 1 is trivial). We will then use an induction to prove the general case. We remark that this type of induction is a standard procedure used in many other places, such as in the context of proving Plücker relations and M ♮ -concavity [MS18].
Before delving into the proof, note that we can always assume S ∩ T = ∅. This is because we can always condition the distribution µ on having any set of elements, and then throwing out those elements; this operation corresponds to taking partial derivatives of g µ which results in a log-concave polynomial by Lemma 11. In particular, we can condition µ on having S ∩ T, and then throwing out S ∩ T from the ground set.
Proof of Lemma 21 for k = 2. When k = 2, we might as well assume that n = 4, because no element outside of S ∪ T is important, and we can condition the distribution µ on not having those elements. This corresponds to substituting 0 for variables outside S ∪ T in g µ which preserves log-concavity.
So our goal now is to show that for a log-concave quadratic polynomial in four variables
we have an O(1)-exchange property. W.l.o.g. assume that S = {1, 2} and T = {3, 4}.
Let us consider ∇ 2 g µ . This is a constant matrix, which has at most one positive eigenvalue by Proposition 15. On the other hand it is a matrix with nonnegative entries, so it must have at least one nonnegative eigenvalue as well. Analyzing the possible signs of the eigenvalues, we see that their product, i.e., the determinant is nonpositive: det(∇ 2 g µ ) ≤ 0.
This determinant can be written in a special way. Let us define:
Notice that approximate exchange for S, T any any i ∈ S is equivalent to saying that A ≤ O(1) • max {B, C}.
We can write det(∇ 2 g µ ) = A 2 + B 2 + C 2 -2(AB + AC + BC). So we get the inequality
Taking square-roots we get
which is the same as saying A ≤ (
Taking square-roots again we get
In particular one of √ B and √ C must be at least 1 2 √ A. This proves that µ satisfies a 2 2 = 4-approximate exchange property for S = {1, 2} and T = {3, 4}.
We now complete the proof by inducting on k.
Proof of Lemma 21 for the general case. We can assume that for any S, T such that |S ∩ T| ≥ 1, we have a 2 O(k-|S∩T|) -approximate exchange property. This is because by the arguments we had, such nonempty intersections can be reduced to smaller values of k by conditioning and throwing out S ∩ T. Now let S ∩ T = ∅ and let i ∈ S be given. Our goal is to find j such that µ(S)µ(T) ≤ 2 O(k) µ(Si + j)µ(T + ij).
Let i ′ ̸ = i be another, arbitrary, element of S. We will exchange i ′ with an element j ′ ∈ T and use induction on Si ′ + j ′ and T. We need to be careful how we choose j ′ though. Let us choose j ′ to be the element of T that maximizes the expression µ(T + ij ′ )µ(Si ′ + j ′ ). The reason for this choice will become apparent in the rest of he proof.
Then the sets Si ′ + j ′ and T have an intersection of one element, so by induction we know an approximate exchange property for them. Therefore, there must be a j ∈ T such that
(4)
We will apply approximate exchange a second time. The sets S and Sii ′ + j + j ′ have a very large intersection. In particular their exchange property reduces to the case of k = 2 of Lemma 21, which we have already proven. By this exchange property, we have
If the first term in Eq. (5) achieves the maximum, then we are done, because multiplying Eqs. (4) and (5) yields
which simplifies to µ(S)µ(T) ≤ 2 O(k) µ(Si + j)µ(T + ij), showing that i can be exchange for j.
So assume that the second term in Eq. (5) achieves the maximum. We will show that in this case i can be exchanged for j ′ . Multiplying Eqs. (4) and (5) yields
Notice that by our choice of j ′ , the fraction appearing on the r.h.s. is ≤ 1. So we can conclude that
If we require the stricter assumption that µ is generated by a real-stable polynomial, then we obtain an exponentially improved k 2 -exchange. This is a generalization of Proposition 23.
Lemma 24. Consider µ : ( [n] k ) → R ≥0 that is generated by a real-stable polynomial. For every S, T ∈ (
Consequently, there exists j ∈ T \ S such that
Thus µ satisfies a k 2 -exchange property. Moreover, for S ∈ ( [n] k ) and j ̸ ∈ S,
where µ(t) = ∑ T∈( [n] k ):t∈T µ(T) for t ∈ {j, e} . A consequence of Lemma 24 is that natural local search and greedy algorithms for finding the maximum of µ : ( [n] k ) → R ≥0 (i.e., MAP inference) give a (k!) 2 ≃ k O(k) -approximation of max {µ(S)}, assuming µ's generating polynomial g µ (z 1 , • • • , z n ) = ∑ S µ(S)z S is real-stable/strongly Rayleigh. This generalizes similar results for determinantal µ to the class of strongly Rayleigh distributions [see, e.g., KD16], giving further evidence for the efficacy of local search and greedy methods [see also Fed13].
Algorithm 2: Greedy Initialize S ← ∅ while |S| < k do Pick i ̸ ∈ S that maximizes µ(S ∪ {i}), and update S ← S ∪ {i} end Lemma 25. If µ is strongly Rayleigh, then the output S ∈ (
k ) µ(T). We leave the proof of Lemmas 24 and 25 to the full version. We remark that similar guarantees can be obtained for a closely related local search algorithm, which moves between sets of size k, each time replacing one element by another. Note that our improved exchange property for strongly Rayleigh distributions is crucial in obtaining k O(k) -approximation. For arbitrary log-concave distributions, we can show the approximate exchange property in Lemma 24 with approximation factor 2 O(k) instead of k 2 , thus proving a 2 O(k 2 ) -approximation guarantee for greedy. Furthermore, we show that local search yields k O(k) -approximation. We leave the proof of these statements to the full version.
In this section we prove Theorems 2 and 3. 3Let µ be the distribution over forests of G defined in Theorem 2. In Lemma 26, we show a homogeneous multiaffine log-concave polynomial f ↑ M,q,w that generates a homogeneous distribution µ ↑ :
and (v ⃗α)e i + e j and (v ⃗β)e j + e i are in supp(g). This implies α j > β j and αe j + e i and βe i + e j are in supp( f ).
We proceed using induction on n. Obviously, for n = 1, f is a linear function in w 0 , w 1 with positive coefficient, which is trivially completely log-concave. Suppose the statement is true for matroids M ′ on ground set [n -1] with n ≥ 2.
We only need to verify ∂ α f is log-concave for all α with |α| = n -2. Note that for i ∈ [n],
is strongly log-concave (and of degree 2, thus real-stable), since it is the sum over the bases of dual matroid M * of M (Theorem 13).
For
where T ⊆ [n] is such that B(M) = {[n] \ i | i ∈ T} are the bases of M, e k is elementary symmetric polynomial of degree k, and z S is shorthand for {z
Note that e 2 is real-stable (Theorem 10), and that u ∈ H whenever (z 0 , z
is the dual or complement of the strongly log-concave polynomial
We remark that while the dual of a real-stable polynomial is real-stable, the dual of a strongly log-concave polynomials is not in general necessarily strongly log-concave.
We are ready to prove Theorem 2.
Proof of Theorem 2. Let M be the graphic matroid on graph G = G(V, E) on n edges. Let k := rank(M). W.l.o.g., we can label the edges by 1, 2 • • • , n and assume E = [n]. Let µ be the distribution over independent sets of M (i.e. forest of G) where µ(F) ∝ q k-|F| w F for F ∈ I(M). Note that we can remove all edge of weight 0 from E without changing µ. W.l.o.g., we assume this is already done; thus w(e) > 0 for all e ∈ E.
Let Y := {y 1 , • • • , y k } , Z := {z 1 , • • • , z n } . We identify the variable z i with the edge labeled by i. Let f := f ↑ M,q,w then f is multi-affine, homogeneous, and completely log-concave by Lemma 26. Observe that f is the generating polynomial for distribution µ ↑ : ( Y∪Z n ) → R ≥0 defined by µ ↑ (T) ∝ q k-|T∩Y| ( k |T∩Y| )
w Z\T if Z \ T ∈ I(M), and 0 otherwise.
We run the down-up walk starting from T 0 ∈ supp(µ ↑ ) (for e.g. T 0 = (Z \ F ) ∪ Y for some spanning forest F ). Let ν ↑ be the distribution of the set T ∈ sup(µ ↑ ) we obtained after O(n log(n/ϵ)) steps; Theorem 1 implies ∥µ ↑ν ↑ ∥ TV ≤ ϵ. We then remove all y j from S i.e. collapse T ⊆ Y ∪ Z to T Z := Z \ (T ∩ Z). Let ν be the distribution of T Z . Clearly, supp(ν) = I(M), and if ν ↑ is the same distribution as µ ↑ , then ν is the same as µ. By the data processing inequality, the total variation between ν and µ is at most ϵ since ∥µ -ν∥ ≤ ∥µ ↑ν ↑ ∥ TV ≤ ϵ.
We show each step of the random walk can be implemented in O(log n) time. In the down step, we keep track of whether the level (n -1) set S still satisfies S Z := Z \ S is a forest. Note that if we dropped an y j in the down step to arrive at S, then S Z is always a forest; if, instead, we dropped a z i (equivalently, added z i to S Z ), then we can check whether S Z stays a forest in O(log n) amortized time using link-cut tree [ST83; RTF18]. If S Z is not a forest, then let cycle S be the unique cycle in S Z which contains the edge z i that was added to S Z . When we perform an up-step from S, if S Z is not a forest, select z f among the edges in cycle S with probability ∝ 1/w f and add it to S (equivalently, remove z f from S Z ). This can be done in O(log n)-amortized time (see [RTF18]). If S Z is empty, then we can only add y j which is not already in S Y := S ∩ Y with uniform probability. If S Z is a nonempty forest, then we can add any variable y j ̸ ∈ S Y or z i ∈ S Z . In this case, the probability of adding variable can be explicitly computed i.e. uniform among y j ̸ ∈ S Y , and P[z i ]/ P[y j ] = q k-(ℓ-1) w -1 i /( k ℓ-1 )
where ℓ := |S Z | = |S Y | + 1. We can perform these operations in O(log n) time by
• With probability 1/(1 + τ) where
sample y j uniformly at random from Y \ S Y and add y j to S Y . Note that this action will always be performed if q = 0.
• Maintain an array of cumulative sums s t := ∑ t h=1 w -1 i h for t ∈ [ℓ] where w i 1 , w i 2 , • • • , w i ℓ are the weights corresponding to the edges in S Z ; this data structure supports amortized O(log n)-time insertion and deletion from S Z and binary search in the sorted array [s t ] ℓ t=1 . This data structure can be implemented using a splay tree where each node stores the sum of all leaves in its rooted subtree. With probability τ/(1 + τ), sample z f from S Z with probability ∝ 1/w f by: sample uniformly random p ∈ [0, s ℓ ], find minimum t ∈ [ℓ] where p ≤ s t , and remove z i t from S Z . This removal will split a tree in the forest S Z , and we update the link-cut tree representation of S Z accordingly in O(log n) time.
See Fig. 1 for a visualization of how one up-then-down step may change the set T Z = Z \ T.
For completeness, we briefly summarize how to handle sampling and removing an edge from cycle S , which was described in [RTF18]. We represent S Z as a forest of link-cut trees. When we add an edge e = (u, v) that forms a cycle, splay u to be the root of its tree T u , then access v (which is also in T u ) so that the entire path P u,v from u to v in T u is stored in one auxiliary tree. This auxiliary splay tree can be augmented to support
Note that the running time is n times the mixing time for the down-up walk.
We note however that one still needs to be able to implement each step of the random walk efficiently when n is infinitely large. For examples where this is possible see[OR18].
We remark that reliance on Theorem 1 in this section is not mandatory and the results of this section would have been possible even without Theorem 1.