{'fraction_non_alphanumeric': 0.04354725502761757, 'fraction_numerical': 0.020171503555189357, 'mean_word_length': 4.2448783610755445, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Zonotope hit-and-run for efficient sampling from projection DPPs

Determinantal point processes (DPPs) are distributions over sets of items that model diversity using kernels. Their applications in machine learning include summary extraction and recommendation systems. Yet, the cost of sampling from a DPP is prohibitive in large-scale applications, which has triggered an effort towards efficient approximate samplers. We build a novel MCMC sampler that combines ideas from combinatorial geometry, linear programming, and Monte Carlo methods to sample from DPPs with a fixed sample cardinality, also called projection DPPs. Our sampler leverages the ability of the hit-and-run MCMC kernel to efficiently move across convex bodies. Previous theoretical results yield a fast mixing time of our chain when targeting a distribution that is close to a projection DPP, but not a DPP in general. Our empirical results demonstrate that this extends to sampling projection DPPs, i.e., our sampler is more sample-efficient than previous approaches which in turn translates to faster convergence when dealing with costly-to-evaluate functions, such as summary extraction in our experiments.


Introduction
Determinantal point processes (DPPs) are distributions over configurations of points that encode diversity through a kernel function. DPPs were introduced by Macchi (1975) and have then found applications in fields as diverse as probability (Hough et al., 2006), number theory (Rudnick & Sarnak, 1996), statistical physics (Pathria & Beale, 2011), Monte Carlo methods (Bardenet & Hardy, 2016), and spatial statistics (Lavancier et al., 2015). In machine learning, DPPs over finite sets have been used as a model of diverse sets of items, where the kernel function takes the form of a finite matrix, see  for a comprehensive survey. Applications of DPPs in machine learning (ML) since this survey also include recommendation tasks (Kathuria et al., 2016;Gartrell et al., 2017), text summarization (Dupuy & Bach, 2016), or models for neural signals (Snoek et al., 2013).
Sampling generic DPPs over finite sets is expensive. Roughly speaking, it is cubic in the number r of items in a DPP sample. Moreover, generic DPPs are sometimes specified through an n × n kernel matrix that needs diagonalizing before sampling, where n is the number of items to pick from. In text summarization, r would be the desired number of sentences for a summary, and n the number of sentences of the corpus to summarize. Thus, sampling quickly becomes intractable for large-scale applications . This has motivated research on fast sampling algorithms. While fast exact algorithms exist for specific DPPs such as uniform spanning trees (Aldous, 1990;Broder, 1989;Propp & Wilson, 1998), generic DPPs have so far been addressed with approximate sampling algorithms, using random projections , low-rank approximations (Kulesza & Taskar, 2011;Gillenwater et al., 2012;Affandi et al., 2013), or using Markov chain Monte Carlo techniques (Kang, 2013;Li et al., 2016a;Rebeschini & Karbasi, 2015;Anari et al., 2016;Li et al., 2016b). In particular, there are polynomial bounds on the mixing rates of natural MCMC chains with arbitrary DPPs as their limiting measure; see Anari et al. (2016) for cardinality-constrained DPPs, and Li et al. (2016b) for the general case.
In this paper, we contribute a non-obvious MCMC chain to approximately sample from projection DPPs, which are DPPs with a fixed sample cardinality. Leveraging a combinatorial geometry result by Dyer & Frieze (1994), we show that sampling from a projection DPP over a finite set can be relaxed into an easier continuous sampling problem with a lot of structure. In particular, the target of this continuous sampling problem is supported on the volume spanned by the columns of the feature matrix associated to the projection DPP, a convex body also called a zonotope. This zonotope can be partitioned into tiles that uniquely correspond to DPP realizations, and the relaxed target distribution is flat on each tile. Previous MCMC approaches to sampling projections DPPs can be viewed as attempting moves between neighboring tiles. Using linear programming, we propose an MCMC chain that moves more freely across this tiling. Our chain is a natural transformation of a fast mixing hit-and-run Markov chain (Lovász & Vempala, 2003) on the underlying zonotope; this empirically results in more uncorrelated MCMC samples than previous work. While the results of Anari et al. (2016) and their generalization by Li et al. (2016b) apply to projection DPPs, our experiments support the fact that our chain mixes faster.
The rest of the paper is organized as follows. In Section 2, we introduce projection DPPs and review existing approaches to sampling. In Section 3, we introduce zonotopes and we tailor the hit-and-run algorithm to our needs. In Section 4, we empirically investigate the performance of our MCMC kernel on synthetic graphs and on a summary extraction task, before concluding in Section 5.

Sampling Projections DPPs
In this section, we introduce projection DPPs in two equivalent ways, respectively following Hough et al. (2006), Lyons (2003). Both definitions shed a different light on the algorithms in Section 3.

Projection DPPs as Particular DPPs
Let E = [n] {1, . . . n}. Let also K be a real symmetric positive semidefinite n × n matrix, and for I ⊂ E, write K I for the square submatrix of K obtained by keeping only rows and columns indexed by I ⊂ E. The random subset X ⊂ E is said to follow a DPP on E = {1, . . . , n} with kernel K if (1) Existence of the DPP described by (1) is guaranteed provided K has all its eigenvalues in [0, 1], see e.g., Kulesza & Taskar (2012, Theorem 2.3). Note that (1) encodes the repulsiveness of DPPs. In particular, for any distinct i, j ∈ [n], In other words, K ij encodes departure from independence. Similarly, for constant K ii , K jj , the larger K 2 ij , the less likely it is to have items i and j co-occur in a sample.
Projection DPPs are the DPPs such that the eigenvalues of K are either 0 or 1, that is, K is the matrix of an orthogonal projection. Projection DPPs are also sometimes called elementary DPPs . One can show that samples from a projection DPP with kernel matrix K almost surely contain r = Tr(K) points and that general DPPs are mixtures of projection DPPs, see e.g., Kulesza & Taskar (2012, Theorem 2.3).

Building Projection DPPs from Linear Matroids
Let r < n, and let A be a full-rank r × n real matrix with columns (a j ) j∈ [n] . . Note that elementary algebra yields that for all B 1 , B 2 ∈ B and x ∈ B 1 \ B 2 , there exists an element y ∈ B 2 \ B 1 such that (3) Property (3) is known as the basis-exchange property. It is used in the definition of general matroids (Oxley, 2003). Lyons (2003) defines a projection DPP as the probability measure on B that assigns to B ∈ B a mass proportional to | det B| 2 , where B A :B is the square matrix formed by the r columns of A indexed by B. Note that this squared determinant is also the squared volume of the parallelotope spanned by the columns indexed by B. In this light, sampling a projection DPP is akin to volume sampling (Deshpande & Rademacher, 2010). Finally, observe that the Cauchy-Binet formula gives the normalization B∈B |det A :B | 2 = det AA T , so that the probability mass assigned to B is gives the equivalence between Sections 2.1 and 2.2.
A fundamental example of DPP defined by a matroid is the random set of edges obtained from a uniform spanning tree (Lyons, 2003). Let G be a connected graph with r + 1 vertices and n edges {e i } i∈ [n] . Let now A be the first r rows of the vertex-edge incidence matrix of G. Then B ⊂ [n] is a basis of M [A] if and only if {e i } i∈B form a spanning tree of G (Oxley, 2003). The transfer current theorem of Burton & Pemantle (1993) implies that the uniform distribution on B is a projection DPP, with kernel matrix (4).

On Projection DPPs and k-DPPs in ML
Projection DPPs are DPPs with realizations of constant cardinality k = r, where r is the rank of K. This constant cardinality is desirable when DPPs are used in summary extraction Dupuy & Bach, 2016) and the size of the required output is predefined. Another way of constraining the cardinality of a DPP is to condition on the event |X| = k, which leads to the so-called k-DPPs . Projection DPPs and k-DPPs are in general different objects. In particular, a k-DPP is not a DPP in the sense of (1) unless its kernel matrix K is a projection. In that sense, k-DPPs are non-DPP objects that generalize projection DPPs. In this paper, we show that projection DPPs can benefit from fast sampling methods. It is not obvious how to generalize our algorithm to k-DPPs.
In ML practice, using projection DPPs is slightly different from using a k-DPP. In some applications, typically with graphs, the DPP is naturally a projection, such as uniform spanning trees described in Section 2.2. But quite often, kernels are built feature-by-feature. That is, for each data item i ∈ [n], a normalized vector of features φ i ∈ R r is chosen, a marginal relevance q i is assigned to item i, and a matrix L is defined as In text summarization, for instance, items i, j could be sentences, q i the marginal relevance of sentence i to the user's query, and φ i features such as tf-idf frequencies of a choice of words, and one could draw from a k-DPP associated to L through P [X = I] ∝ det L I , see e.g., Kulesza & Taskar (2012, Section 4.2.1).
Alternately, let A be the matrix with columns ( and assume r < n and A is full-rank. The latter can be ensured in practice by adding a small i.i.d. Gaussian noise to each entry of A. The projection DPP with kernel K in (4) will yield samples of cardinality r, almost surely, and such that the corresponding columns of A span a large volume, hence feature-based diversity. Thus, if the application requires an output of length p, one can pick r = p, as we do in Section A. Alternatively, if we want an output of size approximately p, we can pick r ≥ p and independently thin the resulting sample, which preserves the DPP structure (Lavancier et al., 2015).

Exact Sampling of Projection DPPs
Hough et al. (2006) give an algorithm to sample general DPPs, which is based on a subroutine to sample projection DPPs. Consider a projection DPP with kernel K such that Tr(K) = r, Hough et al.'s (2006) algorithm follows the chain rule to sample a vector (x 1 , . . . , x r ) ∈ [n] r with successive conditional densities  (4), through QR decomposition of A T . Still, depending on n and r, sampling can become intractable. This has sparked interest in fast approximate sampling methods for DPPs, which we survey in Section 2.5.
Interestingly, there exist fast and exact methods for sampling some specific DPPs, which are not based on the approach of Hough et al. (2006). We introduced the DPP behind uniform spanning trees on a connected graph G in Section 2.2. Random walk algorithms such as the ones by Aldous (1990), Broder (1989), and Propp & Wilson (1998) sample uniform spanning trees in time bounded by the cover time of the graph, for instance, which is O(r 3 ) and can be o(r 3 ) (Levin et al., 2009), where G has r + 1 vertices. This compares favorably with the algorithm of Hough et al. (2006) above, since each sample contains r edges. The Aldous-Broder algorithm, for instance, starts from an empty set T = ∅ and an arbitrary node x 0 , and samples a simple random walk (X t ) t∈N on the edges of G, starting from X 0 = x 0 , and adding edge [X t , X t+1 ] to T the first time it visits vertex X t+1 . The algorithm stops when each vertex has been seen at least once, that is, at the cover time of the graph.

Approximate Sampling of Projection DPPs
There are two main sets of methods for approximate sampling from general DPPs. The first set uses the generalpurpose tools from numerical algebra and the other is based on MCMC sampling.
Consider K = C T C with C of size d × n, for some d n (Kulesza & Taskar, 2011), but still too large for exact sampling using the method of Hough et al. (2006), then Gillenwater et al. (2012) show how projecting C can give an approximation with bounded error. When this decomposition of the kernel is not possible, Affandi et al. (2013) adapt Nyström sampling (Williams & Seeger, 2001) to DPPs and bound the approximation error for DPPs and k-DPPs, which thus applies to projection DPPs.
Apart from general purpose approximate solvers, there exist MCMC-based methods for approximate sampling from projection DPPs. In Section 2.2, we introduced the basisexchange property, which implies that once we remove an element from a basis B 1 of a linear matroid, any other basis B 2 has an element we can take and add to B 1 to make it a basis again. This means we can construct a connected Algorithm 1 basisExchangeSampler Input: Either A or K Initialize i ← 0 and pick B 0 ∈ B as defined in (2) while Not converged do graph G be with B as vertex set, and we add an edge between two bases if their symmetric difference has cardinality 2. G be is called the basis-exchange graph. Feder & Mihail (1992) show that the simple random walk on G be has limiting distribution the uniform distribution on B and mixes fast, under conditions that are satisfied by the matroids involved by DPPs.
If the uniform distribution on B is not the DPP we want to sample from, 1 we can add an accept-reject step after each move to make the desired DPP the limiting distribution of the walk. Adding such an acceptance step and a probability to stay at the current basis, Anari et al. (2016); Li et al. (2016b) give precise polynomial bounds on the mixing time of the resulting Markov chains. This Markov kernel on B is given in Algorithm 1. Note that we use the acceptance ratio of Li et al. (2016b). In the following, we make use of the notation Vol defined as follows. For any P ⊂ [n], which corresponds to the squared volume of the parallelotope spanned by the columns of A indexed by P . In particular, for subsets P such that |P | > r or such that |P | = r, P / ∈ B we have Vol 2 (A :P ) = 0. However, for B ∈ B, We now turn to our contribution, which finds its place in this category of MCMC-based approximate DPP samplers.

Hit-and-run on Zonotopes
Our main contribution is the construction of a fast-mixing Markov chain with limiting distribution a given projection DPP. Importantly, we assume to know A in (4).
Assumption 1. We know a full-rank r × n matrix A such that K = A T (AA T ) −1 A.
As discussed in Section 2.3, this is not an overly restrictive assumption, as many ML applications start with building the feature matrix A rather than the similarity matrix K.

Zonotopes
We define the zonotope Z(A) of A as the r-dimensional volume spanned by the column vectors of A, As an affine transformation of the unit hypercube, Z(A) is a r-dimensional polytope. In particular, for a basis B ∈ B of the matroid M [A], the corresponding Z(B) is a rdimensional parallelotope with volume Vol(B) = |det B|, see Figure 1(a). On the contrary, any P ⊂ [n], such that |P | = r, P / ∈ B also yields a parallelotope Z(A :P ), but its volume is null. In the latter case, the exchange move in Algorithm 1 will never be accepted and the state space of the corresponding Markov chain is indeed B.
Our algorithm relies on the proof of the following.

Vol(Z(
Proof. In short, for a good choice of c ∈ R n , Dyer & Frieze (1994) consider for any x ∈ Z(A), the following linear program (LP) noted P x (A, c), Standard LP results (Luenberger & Ye, 2008) yield that the unique optimal solution y * of P x (A, c) takes the form with u ∈ [0, 1] r and ξ(x) ∈ {0, 1} n such that ξ(x) i = 0 for i ∈ B x . In case the choice of B x is ambiguous, Dyer & Frieze (1994) take the smallest in the lexicographic order. Decomposition (10) allows locating any point x ∈ Z(A) as falling inside a uniquely defined parallelotope Z(B x ) shifted by ξ(x). Manipulating the optimality conditions of (9), Dyer & Frieze (1994) prove that each basis B can be realized as a B x for some x, and that  Note that c is used to fix the tiling of the zonotope, but the map x → B x depends on this linear objective. Therefore, the tiling of Z(A) is may not be unique. An arbitrary c gives a valid tiling, as long as there are no ties when solving (9). Dyer & Frieze (1994) use a nonlinear mathematical trick to fix c. In practice (Section 4.1), we generate a random Gaussian c once and for all, which makes sure no ties appear during the execution, with probability 1. Remark 1. We propose to interpret the proof of Proposition 1 as a volume sampling algorithm: if one manages to sample an x uniformly on Z(A), and then extracts the corresponding basis B = B x by solving (9), then B is drawn with probability proportional to Vol(B) = | det B|.

Remark 1 is close to what we want, as sampling from a projection DPP under Assumption 1 boils down to sampling a basis B of M [A]
proportionally to the squared volume | det B| 2 (Section 2.2). In the rest of this section, we explain how to efficiently sample x uniformly on Z(A), and how to change the volume into its square.

Hit-and-run and the Simplex Algorithm
Z(A) is a convex set. Approximate uniform sampling on large-dimensional convex bodies is one of the core questions in MCMC, see e.g., Cousins & Vempala (2016) and references therein. The hit-and-run Markov chain (Turčin, 1971;Smith, 1984) is one of the preferred practical and theoretical solutions (Cousins & Vempala, 2016).
We describe the Markov kernel P (x, z) of the hit-and-run Markov chain for a generic target distribution π supported on a convex set C. Sample a point y uniformly on the unit sphere centered at x. Letting d = y − x, this defines the line D x {x + αd ; α ∈ R}. Then, sample z from any Markov kernel Q(x, ·) supported on D x that leaves the restriction of π to D x invariant. In particular, Metropolis-Hastings kernel (MH, Robert & Casella 2004) is often used with uniform proposal on D x , which favors large moves across the support C of the target, see Figure 1(b). The resulting Markov kernel leaves π invariant, see e.g., Andersen & Diaconis (2007) for a general proof. Furthermore, the hit-and-run Markov chain has polynomial mixing time for log concave π (Lovász & Vempala, 2003, Theorem 2.1).
To implement Remark 1, we need to sample from π u ∝ 1 Z(A) . In practice, we can choose the secondary Markov kernel Q(x, ·) to be MH with uniform proposal on D x , as long as we can determine the endpoints x + α m (y − x) and x + α M (y − x) of D x ∩ Z(A). In fact, zonotopes are tricky convex sets, as even an oracle saying whether a point belongs to the zonotope requires solving LPs (basically, it is Phase I of the simplex algorithm). As noted by Lovász & Vempala (2003, Section 4.4), hit-and-run with LP is the state-of-the-art for computing the volume of large-scale zonotopes. Thus, by definition of Z(A), this amounts to solving two more LPs: α m is the optimal solution to the linear program min λ∈R n ,α∈R α s.t.
x + αd = Aλ 0 ≤ λ ≤ 1, while α M is the optimal solution of the same linear program with objective −α. Thus, a combination of hit-andrun and LP solvers such as Dantzig's simplex algorithm (Luenberger & Ye, 2008) yields a Markov kernel with invariant distribution 1 Z(A) , summarized in Algorithm 2.
Algorithm 2 unifZonoHitAndRun Input: A Initialization: The acceptance in MH is 1 due to our choice of the proposal and the target. By the proof of Proposition 1, running Algorithm 2, taking the output chain (x i ) and extracting the bases (B xi ) with Algorithm 3, we obtain a chain on B with invariant distribution proportional to the volume of B.
In terms of theoretical performance, this Markov chain inherits Lovász & Vempala's (2003) mixing time as it is a simple transformation of hit-and-run targeting the uniform distribution on a convex set. We underline that this is not a pathological case and it already covers a range of applications, as changing the feature matrix A yields another zonotope, but the target distribution on the zonotope stays uniform. Machine learning practitioners do not use volume sampling for diversity sampling yet, but nothing prevents it, as it already encodes the same feature-based diversity as squared volume sampling (i.e., DPPs). Nevertheless, our initial goal was to sample from a projection DPP with kernel K under Assumption 1. We now modify the Markov chain just constructed to achieve that.

From Volume to Squared Volume
Consider the probability density function on Z(A) represented on our example in Figure 1(c). Observe, in particular, that π v is constant on each Z(B). Running the hit-and-run algorithm with this target instead of π u in Section 3.2, and extracting bases using Algorithm 3 again, we obtain a Markov chain on B with limiting distribution ν(B) proportional to the squared volume spanned by column vectors of B, as required. To see this, note that ν(B) is the volume of the "skyscraper" built on top of Z(B) in Figure 1(c), that is Vol(B) × Vol(B).

Algorithm 4 volZonoHitAndRun
The resulting algorithm is shown in Algorithm 5. Note the acceptance ratio in the subroutine Algorithm 4 compared to Algorithm 2, since the target of the hit-and-run algorithm is not uniform anymore.

On Base Measures
As described in Section 2.3, it is common in ML to specify a marginal relevance q i of each item i ∈ [n], i.e., the base measure of the DPP. Compared to a uniform base measure, this means replacing A by A with columns a i = √ q i a i .
Contrary to A, in Algorithm 4, both the zonotope and the acceptance ratio are scaled by the corresponding products of √ q i s. We could equally well define A by multiplying each column of A by q i instead of its square root, and leave the acceptance ratio in Algorithm 4 use columns of the original A. By the arguments in Section 3.3, the chain (B i ) would leave the same projection DPP invariant. In particular, we have some freedom in how to introduce the marginal relevance q i , so we can choose the latter solution that simply scales the zonotope and its tiles to preserve outer angles, while using unscaled volumes to decide acceptance. This way, we do not create harder-to-escape or sharper corners for hit-and-run, which could lead the algorithm to be stuck for a while (Cousins & Vempala, 2016, Section 4.2.1). Finally, since hit-and-run is efficient at moving across convex bodies (Lovász & Vempala, 2003), the rationale is that if hit-and-run was empirically mixing fast before scaling, its performance should not decrease.

Experiments
We investigate the behavior of our Algorithm 5 on synthetic graphs in Section 4.1, in summary extraction in Section 4.2, and on MNIST in Appendix A.

Non-uniform Spanning Trees
We compare Algorithm 1 studied by Anari et al. (2016); Li et al. (2016b) and our Algorithm 5 on two types of graphs, in two different settings. The graphs we consider are the complete graph K 10 with 10 vertices (and 45 edges) and a realization BA(20, 2) of a Barabási-Albert graph with 20 vertices and parameter 2. We chose BA as an example of structured graph, as it has the preferential attachment property present in social networks (Barabási & Albert, 1999). The input matrix A is a weighted version of the vertex-edge incidence matrix of each graph for which we keep only the 9 (resp. 19) first rows, so that it satisfies Assumption 1. For more generality, we introduce a base measure, as described in Section 2.3 and 3.4, by reweighting the columns of A with i.i.d. uniform variables in [0, 1]. Samples from the corresponding projection DPP are thus spanning trees drawn proportionally to the products of their edge weights.
For Algorithm 5, a value of the linear objective c is drawn once and for all, for each graph, from a standard Gaussian distribution. This is enough to make sure no ties appear during the execution, as mentioned in Section 3.1. This linear objective is kept fixed throughout the experiments so that the tiling of the zonotope remains the same. We run both algorithms for 70 seconds, which corresponds to roughly 50 000 iterations of Algorithm 5. Moreover, we run 100 chains in parallel for each of the two algorithms. For each of the 100 repetitions, we initialize the two algorithms with the same random initial basis, obtained by solving (9) once, with x = Au and u ∼ U [0,1] n . For both graphs, the total number |B| of bases is of order 10 8 , so computing total variation distances is impractical. We instead compare Algorithms 1 and 5 based on the estimation of inclusion probabilities P [S ⊂ B] for various subsets S ⊂ [n] of size 3. We observed similar behaviors across 3-subsets, so we display here the typical behavior on a 3-subset.
The inclusion probabilities are estimated via a running average of the number of bases containing the subsets S. show the behavior of both algorithms vs. wall-clock time for the complete graph K 10 and a realization of BA(20, 2), respectively. In these four figures, bold curves correspond to the median of the relative errors, whereas the frontiers of colored regions indicate the first and last deciles of the relative errors.
In Figures 2(c) and 3(c) we compute the Gelman-Rubin statistic (Gelman & Rubin, 1992), also called the potential scale reduction factor (PSRF). We use the PSRF implementation of CODA (Plummer et al., 2006) in R, on the 100 binary chains indicating the presence of the typical 3-subset in the current basis.
In terms of number of iterations, our Algorithm 5 clearly mixes faster. Relatedly, we observed typical acceptance rates for our algorithm an order of magnitude larger than Algorithm 1, while simultaneously attempting more global moves than the local basis-exchange moves of Algorithm 1. The high acceptance is partly due to the structure of the zonotope: the uniform proposal in the hit-and-run algorithm already favors bases with large determinants, as the length of the intersection of D x in Algorithm 4 with any Z(B) is an indicator of its volume, see also Figure 1 Under the time-horizon constraint, see Figures 2(b) and 3(b), Algorithm 1 has time to perform more than 10 6 iterations compared to roughly 50 000 steps for our chain. The acceptance rate of Algorithm 5 is still 10 times larger, but the time required to solve the linear programs at each MCMC iteration clearly hinders our algorithm in terms of CPU time. Both algorithms are comparable in performance, but given its large acceptance, we would expect our algorithm to perform better if it was allowed to do even only 10 times more iterations. Now this is implementationdependent, and our current implementation of Algorithm 5 is relatively naive, calling the simplex algorithm in the GLPK (Oki, 2012) solver with CVXOPT (Andersen et al., 2008) from Python. We think there are big potential speedups to realize in the integration of linear programming solvers in our code. Moreover, we initialize our simplex algorithms randomly, while the different LPs we solve are related, so there may be additional smart mathematical speed-ups in using the path followed by one simplex instance to initialize the next.
Finally, we note that the performance of our Algorithm 5 seems stable and independent of the structure of the graph, while the performance of the basis-exchange Algorithm 1 seems more graph-dependent. Further investigation is needed to make stronger statements.

Text Summarization
Looking at Figures 2 and 3, our algorithm will be most useful when the bottleneck is mixing vs. number of iterations rather than CPU time. For instance, when integrating a costly-to-evaluate function against a projection DPP, the evaluation of the integrand may outweigh the cost of one iteration. To illustrate this, we adapt an experiment of Kulesza & Taskar (2012, Section 4.2.1) on minimum Bayes risk decoding for summary extraction. The idea is to find a   subset Y of sentences of a text that maximizes where (Y r ) r are sampled from a projection DPP. ROUGE-1F is a measure of similarity of two sets of sentences. We summarize this 64-sentence article as a subset of 11 sentences. In this setting, evaluating once ROUGE-1F in the sum (12) takes 0.1s on a modern laptop, while one iteration of our algorithm is 10 −3 s. Our Algorithm 5 can thus compute (12) for R = 10 000 in about the same CPU time as Algorithm 1, an iteration of which costs 10 −5 s. We show in Figure 4 the value of (12) for 3 possible summaries Y (i) 3 i=1 chosen uniformly at random in B, over 50 independent runs. The variance of our estimates is smaller, and the number of different summaries explored is about 50%, against 10% for Algorithm 1. Evaluating (12) using our algorithm is thus expected to be closer to the maximum of the underlying integral. Details are given in Appendix B.

Discussion
We proposed a new MCMC kernel with limiting distribution being an arbitrary projection DPP. This MCMC kernel leverages optimization algorithms to help making global moves on a convex body that represents the DPP. We provided empirical results supporting its fast mixing when compared to the state-of-the-art basis-exchange chain of Anari et al. (2016); Li et al. (2016b). Future work will focus on an implementation: while our MCMC chain mixes faster, when compared based on CPU time our algorithm suffers from having to solve linear programs at each iter- ation. We note that even answering the question whether a given point belongs to a zonotope involves linear programming, so that chord-finding procedures used in slice sampling (Neal, 2003, Sections 4 and 5) would not provide significant computational savings.
We also plan to investigate theoretical bounds on the mixing time of our Algorithm 4. We can build upon the work of Anari et al. (2016), as our Algorithm 4 is also a weighted extension of our Algorithm 2, and the polynomial bounds for the vanilla hit-and-run algorithm (Lovász & Vempala, 2003) already apply to the latter. Note that while not targeting a DPP, our Algorithm 2 already samples items with feature-based repulsion, and could be used independently if the determinantal aspect is not crucial to the application.

A. Landmarks on MNIST
To illustrate our Algorithm 5 on a non-synthetic dataset, we take a thousand images of handwritten 1s from MNIST. We consider the problem of finding landmarks among these 1s, see, for example, Liang & Paisley (2015). To obtain a set of r = 10 landmarks from a projection DPP, we design r-dimensional feature vectors of our images, the inner products of which indicate similarity. Specifically, we take A to be the coordinates of the projections of our dataset onto the first r principal directions given by a standard PCA, so that A has rank r if the whole variance is not explained by these first r eigenvalues. Note that our choice of using standard PCA is arbitrary: For manifold landmarking, we could take the coordinates output by a manifold learning algorithm such as ISOMAP.
Running our Algorithm 5 for 10 000 time steps, we obtain bases with squared volumes spanning three orders of magnitude. The basis in our MCMC sample with maximum squared volume is shown in the first row of Figure 5. The bottom three rows show bases drawn uniformly at random from our initial dataset, for visual comparison. The left column gives the log of the ratio of the squared volume of the corresponding basis by that of the top row.  Visual inspection reveals that the projection DPP and our MCMC sample successfully pick up diversity among 1s: Different angles, thicknesses, and shapes are selected, while uniform random sampling exhibits less coverage. This is confirmed by the log ratios, which are far from the three orders of magnitude within our chain.

B. Extractive Text Summarization
We consider the article 2 entitled Scientists, Stop Thinking Explaining Science Will Fix Things. In order to generate an 11-sentences summary, we build 11 features as follows. For each sentence, we compute its number of characters and its number of words. Then, we apply a Porter stemmer (Loper & Bird, 2002) and count again the number of characters and words in each sentence. In addition, we sum the tf-idf values of the words in each sentence and compute the average cosine distance to all other sentences. Finally, we compute the position of the sentence in the original article and generate binary features indicating positions 15. We end up with a feature matrix A of size 11 × 64.
Next, we display 3 summaries of the article, the first one is drawn using Algorithm 5 while the 2 others are constructed by picking 11 sentences uniformly at random. In both cases, we sort the sentences in the order they originally appear. Samples from Algorithm 5 induce a better coverage of the article than uniform draws.