Determinantal point processes for machine learning

Introduction

Probabilistic modeling and learning techniques have become indispensable tools for analyzing data, discovering patterns, and making predictions in a variety of real-world settings. In recent years, the widespread availability of both data and processing capacity has led to new applications and methods involving more complex, structured output spaces, where the goal is to simultaneously make a large number of interrelated decisions. Unfortunately, the introduction of structure typically involves a combinatorial explosion of output possibilities, making inference computationally impractical without further assumptions.

A popular compromise is to employ graphical models, which are tractable when the graph encoding local interactions between variables is a tree. For loopy graphs, inference can often be approximated in the special case when the interactions between variables are positive and neighboring nodes tend to have the same labels. However, dealing with global, negative interactions in graphical models remains intractable, and heuristic methods often fail in practice.

Determinantal point processes (DPPs) offer a promising and complementary approach. Arising in quantum physics and random matrix theory, DPPs are elegant probabilistic models of global, negative correlations, and offer efficient algorithms for sampling, marginalization, conditioning, and other inference tasks. While they have been studied extensively by mathematicians, giving rise to a deep and beautiful theory, DPPs are relatively new in machine learning. We aim to provide a comprehensible introduction to DPPs, focusing on the intuitions, algorithms, and extensions that are most relevant to our community.

A DPP is a distribution over subsets of a fixed ground set, for instance, sets of search results selected from a large database. Equivalently, a DPP over a ground set of $N$ items can be seen as modeling a binary characteristic vector of length $N$. The essential characteristic of a DPP is that these binary variables are negatively correlated; that is, the inclusion of one item makes the inclusion of other items less likely. The strengths of these negative correlations are derived from a kernel matrix that defines a global measure of similarity between pairs of items, so that more similar items are less likely to co-occur. As a result, DPPs assign higher probability to sets of items that are diverse; for example, a DPP will prefer search results that cover multiple distinct aspects of a user’s query, rather than focusing on the most popular or salient one.

This focus on diversity places DPPs alongside a number of recently developed techniques for working with diverse sets, particularly in the information retrieval community (Carbonell and Goldstein, 1998; Zhai et al., 2003; Chen and Karger, 2006; Yue and Joachims, 2008; Radlinski et al., 2008; Swaminathan et al., 2009; Raman et al., 2012). However, unlike these methods, DPPs are fully probabilistic, opening the door to a wider variety of potential applications, without compromising algorithmic tractability.

The general concept of diversity can take on a number of forms depending on context and application. Including multiple kinds of search results might be seen as covering or summarizing relevant interpretations of the query or its associated topics; see Figure 1. Alternatively, items inhabiting a continuous space may exhibit diversity as a result of repulsion, as in Figure 2. In fact, certain repulsive quantum particles are known to be precisely described by a DPP; however, a DPP can also serve as a model for general repulsive phenomena, such as the locations of trees in a forest, which appear diverse due to physical and resource constraints. Finally, diversity can be used as a filtering prior when multiple selections must be based on a single detector or scoring metric. For instance, in Figure 3 a weak pose detector favors large clusters of poses that are nearly identical, but filtering through a DPP ensures that the final predictions are well-separated.

Throughout this survey we demonstrate applications for DPPs in a variety of settings, including:

The DUC 2003/2004 text summarization task, where we form extractive summaries of news articles by choosing diverse subsets of sentences (Section 4.2.1);

An image search task, where we model human judgments of diversity for image sets returned by Google Image Search (Section 5.3);

A multiple pose estimation task, where we improve the detection of human poses in images from television shows by incorporating a bias toward non-overlapping predictions (Section 6.4);

A news threading task, where we automatically extract timelines of important news stories from a large corpus by balancing intra-timeline coherence with inter-timeline diversity (Section 6.6.4).

2 Outline

In this paper we present general mathematical background on DPPs along with a range of modeling extensions, efficient algorithms, and theoretical results that aim to enable practical modeling and learning. The material is organized as follows.

We begin with an introduction to determinantal point processes tailored to the interests of the machine learning community. We focus on discrete DPPs, emphasizing intuitions and including new, simplified proofs for some theoretical results. We provide descriptions of known efficient inference algorithms, and characterize their computational properties.

We describe a decomposition of the DPP that makes explicit its fundamental tradeoff between quality and diversity. We compare the expressive power of DPPs and MRFs, characterizing the tradeoffs in terms of modeling power and computational efficiency. We also introduce a dual representation for DPPs, showing how it can be used to perform efficient inference over large ground sets. When the data are high-dimensional and dual inference is still too slow, we show that random projections can be used to maintain a provably close approximation to the original model while greatly reducing computational requirements.

We derive an efficient algorithm for learning the parameters of a quality model when the diversity model is held fixed. We employ this learning algorithm to perform extractive summarization of news text.

We present an extension of DPPs that allows for explicit control over the number of items selected by the model. We show not only that this extension solves an important practical problem, but also that it increases expressive power: a $k$-DPP can capture distributions that a standard DPP cannot. The extension to $k$-DPPs necessitates new algorithms for efficient inference based on recursions for the elementary symmetric polynomials. We validate the new model experimentally on an image search task.

We extend DPPs to model diverse sets of structured items, such as sequences or trees, where there are combinatorially many possible configurations. In this setting the number of possible subsets is doubly-exponential, presenting a daunting computational challenge. However, we show that a factorization of the quality and diversity models together with the dual representation for DPPs makes efficient inference possible using second-order message passing. We demonstrate structured DPPs on a toy geographical paths problem, a still-image multiple pose estimation task, and two high-dimensional text threading tasks.

Determinantal point processes

Determinantal point processes (DPPs) were first identified as a class by Macchi (1975), who called them “fermion processes” because they give the distributions of fermion systems at thermal equilibrium. The Pauli exclusion principle states that no two fermions can occupy the same quantum state; as a consequence fermions exhibit what is known as the “anti-bunching” effect. This repulsion is described precisely by a DPP.

In fact, years before Macchi gave them a general treatment, specific DPPs appeared in major results in random matrix theory (Mehta and Gaudin, 1960; Dyson, 1962a, b, c; Ginibre, 1965), where they continue to play an important role (Diaconis, 2003; Johansson, 2005b). Recently, DPPs have attracted a flurry of attention in the mathematics community (Borodin and Olshanski, 2000; Borodin and Soshnikov, 2003; Borodin and Rains, 2005; Borodin et al., 2010; Burton and Pemantle, 1993; Johansson, 2002, 2004, 2005a; Okounkov, 2001; Okounkov and Reshetikhin, 2003; Shirai and Takahashi, 2000), and much progress has been made in understanding their formal combinatorial and probabilistic properties. The term “determinantal” was first used by Borodin and Olshanski (2000), and has since become accepted as standard. Many good mathematical surveys are now available (Borodin, 2009; Hough et al., 2006; Shirai and Takahashi, 2003a, b; Lyons, 2003; Soshnikov, 2000; Tao, 2009).

We begin with an overview of the aspects of DPPs most relevant to the machine learning community, emphasizing intuitions, algorithms, and computational properties.

A point process $\mathcal{P}$ on a ground set $\mathcal{Y}$ is a probability measure over “point patterns” or “point configurations” of $\mathcal{Y}$, which are finite subsets of $\mathcal{Y}$. For instance, $\mathcal{Y}$ could be a continuous time interval during which a scientist records the output of a brain electrode, with $\mathcal{P}(\{y_{1},y_{2},y_{3}\})$ characterizing the likelihood of seeing neural spikes at times $y_{1}$, $y_{2}$, and $y_{3}$. Depending on the experiment, the spikes might tend to cluster together, or they might occur independently, or they might tend to spread out in time. $\mathcal{P}$ captures these correlations.

For the remainder of this paper, we will focus on discrete, finite point processes, where we assume without loss of generality that $\mathcal{Y}=\{1,2,\dots,N\}$; in this setting we sometimes refer to elements of $\mathcal{Y}$ as items. Much of our discussion extends to the continuous case, but the discrete setting is computationally simpler and often more appropriate for real-world data—e.g., in practice, the electrode voltage will only be sampled at discrete intervals. The distinction will become even more apparent when we apply our methods to $\mathcal{Y}$ with no natural continuous interpretation, such as the set of documents in a corpus.

In the discrete case, a point process is simply a probability measure on $2^{\mathcal{Y}}$, the set of all subsets of $\mathcal{Y}$. A sample from $\mathcal{P}$ might be the empty set, the entirety of $\mathcal{Y}$, or anything in between. $\mathcal{P}$ is called a determinantal point process if, when $\boldsymbol{Y}$ is a random subset drawn according to $\mathcal{P}$, we have, for every $A\subseteq\mathcal{Y}$,

for some real, symmetric $N\times N$ matrix $K$ indexed by the elements of $\mathcal{Y}$.In general, $K$ need not be symmetric. However, in the interest of simplicity, we proceed with this assumption; it is not a significant limitation for our purposes. Here, $K_{A}\equiv\left[K_{ij}\right]_{i,j\in A}$ denotes the restriction of $K$ to the entries indexed by elements of $A$, and we adopt $\det(K_{\emptyset})=1$. Note that normalization is unnecessary here, since we are defining marginal probabilities that need not sum to 1.

Since $\mathcal{P}$ is a probability measure, all principal minors $\det(K_{A})$ of $K$ must be nonnegative, and thus $K$ itself must be positive semidefinite. It is possible to show in the same way that the eigenvalues of $K$ are bounded above by one using Equation (35), which we introduce later. These requirements turn out to be sufficient: any $K$, $0\preceq K\preceq I$, defines a DPP. This will be a consequence of Theorem 2.3.

We refer to $K$ as the marginal kernel since it contains all the information needed to compute the probability of any subset $A$ being included in $\boldsymbol{Y}$. A few simple observations follow from Equation (1). If $A=\{i\}$ is a singleton, then we have

That is, the diagonal of $K$ gives the marginal probabilities of inclusion for individual elements of $\mathcal{Y}$. Diagonal entries close to 1 correspond to elements of $\mathcal{Y}$ that are almost always selected by the DPP. Furthermore, if $A=\{i,j\}$ is a two-element set, then

Thus, the off-diagonal elements determine the negative correlations between pairs of elements: large values of $K_{ij}$ imply that $i$ and $j$ tend not to co-occur.

Equation (7) demonstrates why DPPs are “diversifying”. If we think of the entries of the marginal kernel as measurements of similarity between pairs of elements in $\mathcal{Y}$, then highly similar elements are unlikely to appear together. If $K_{ij}=\sqrt{K_{ii}K_{jj}}$, then $i$ and $j$ are “perfectly similar” and do not appear together almost surely. Conversely, when $K$ is diagonal there are no correlations and the elements appear independently. Note that DPPs cannot represent distributions where elements are more likely to co-occur than if they were independent: correlations are always nonpositive.

Figure 4 shows the difference between sampling a set of points in the plane using a DPP (with $K_{ij}$ inversely related to the distance between points $i$ and $j$), which leads to a relatively uniformly spread set with good coverage, and sampling points independently, which results in random clumping.

In this paper, our focus is on using DPPs to model real-world data. However, many theoretical point processes turn out to be exactly determinantal, which is one of the main reasons they have received so much recent attention. In this section we briefly describe a few examples; some of them are quite remarkable on their own, and as a whole they offer some intuition about the types of distributions that are realizable by DPPs. Technical details for each example can be found in the accompanying reference.

Given a sequence of $N$ random numbers drawn uniformly and independently from a finite set (say, the digits 0–9), the locations in the sequence where the current number is less than the previous number form a subset of $\{2,3,\dots,N\}$. This subset is distributed as a determinantal point process. Intuitively, if the current number is less than the previous number, it is probably not too large, thus it becomes less likely that the next number will be smaller yet. In this sense, the positions of decreases repel one another.

Let $G$ be an arbitrary finite graph with $N$ edges, and let $T$ be a random spanning tree chosen uniformly from the set of all the spanning trees of $G$. The edges in $T$ form a subset of the edges of $G$ that is distributed as a DPP. The marginal kernel in this case is the transfer-impedance matrix, whose entry $K_{e_{1}e_{2}}$ is the expected signed number of traversals of edge $e_{2}$ when a random walk begins at one endpoint of $e_{1}$ and ends at the other (the graph edges are first oriented arbitrarily). Thus, edges that are in some sense “nearby” in $G$ are similar according to $K$, and as a result less likely to participate in a single uniformly chosen spanning tree. As this example demonstrates, some DPPs assign zero probability to sets whose cardinality is not equal to a particular $k$; in this case, $k$ is the number of nodes in the graph minus one—the number of edges in any spanning tree. We will return to this issue in Section 5.

Let $M$ be a random matrix obtained by drawing every entry independently from the complex normal distribution. This is the complex Ginibre ensemble. The eigenvalues of $M$, which form a finite subset of the complex plane, are distributed according to a DPP. If a Hermitian matrix is generated in the corresponding way, drawing each diagonal entry from the normal distribution and each pair of off-diagonal entries from the complex normal distribution, then we obtain the Gaussian unitary ensemble, and the eigenvalues are now a DPP-distributed subset of the real line.

The Aztec diamond is a diamond-shaped union of lattice squares, as depicted in Figure 5a. (Half of the squares have been colored gray in a checkerboard pattern.) A domino tiling is a perfect cover of the Aztec diamond using $2\times 1$ rectangles, as in Figure 5b. Suppose that we draw a tiling uniformly at random from among all possible tilings. (The number of tilings is known to be exponential in the width of the diamond.) We can identify this tiling with the subset of the squares that are (a) painted gray in the checkerboard pattern and (b) covered by the left half of a horizontal tile or the bottom half of a vertical tile (see Figure 5c). This subset is distributed as a DPP.

2 L-ensembles

For the purposes of modeling real data, it is useful to slightly restrict the class of DPPs by focusing on L-ensembles. First introduced by Borodin and Rains (2005), an L-ensemble defines a DPP not through the marginal kernel $K$, but through a real, symmetric matrix $L$ indexed by the elements of $\mathcal{Y}$:

Whereas Equation (1) gave the marginal probabilities of inclusion for subsets $A$, Equation (8) directly specifies the atomic probabilities for every possible instantiation of $\boldsymbol{Y}$. As for $K$, it is easy to see that $L$ must be positive semidefinite. However, since Equation (8) is only a statement of proportionality, the eigenvalues of $L$ need not be less than one; any positive semidefinite $L$ defines an L-ensemble. The required normalization constant can be given in closed form due to the fact that $\sum_{Y\subseteq\mathcal{Y}}\det(L_{Y})=\det(L+I)$, where $I$ is the $N\times N$ identity matrix. This is a special case of the following theorem.

where $I_{\bar{A}}$ is the diagonal matrix with ones in the diagonal positions corresponding to elements of $\bar{A}=\mathcal{Y}-A$, and zeros everywhere else.

Suppose that $A=\mathcal{Y}$; then Equation (9) holds trivially. Now suppose inductively that the theorem holds whenever $\bar{A}$ has cardinality less than $k$. Given $A$ such that $|\bar{A}|=k>0$, let $i$ be an element of $\mathcal{Y}$ where $i\in\bar{A}$. Splitting blockwise according to the partition $\mathcal{Y}=\{i\}\cup\mathcal{Y}-\{i\}$, we can write

where $L_{\bar{i}i}$ is the subcolumn of the $i$th column of $L$ whose rows correspond to $\bar{i}$, and similarly for $L_{i\bar{i}}$. By multilinearity of the determinant, then,

We can now apply the inductive hypothesis separately to each term, giving

where we observe that every $Y$ either contains $i$ and is included only in the first sum, or does not contain $i$ and is included only in the second sum. ∎

As a shorthand, we will write $\mathcal{P}_{L}(Y)$ instead of $\mathcal{P}_{L}(\boldsymbol{Y}=Y)$ when the meaning is clear.

We can write a version of Equation (7) for L-ensembles, showing that if $L$ is a measure of similarity then diversity is preferred:

In this case we are reasoning about the full contents of $\boldsymbol{Y}$ rather than its marginals, but the intuition is essentially the same. Furthermore, we have the following result of Macchi (1975).

An L-ensemble is a DPP, and its marginal kernel is

Using Theorem 2.1, the marginal probability of a set $A$ is

We can use the fact that $L(L+I)^{-1}=I-(L+I)^{-1}$ to simplify and obtain

where we let $K=I-(L+I)^{-1}$. Now, we observe that left multiplication by $I_{A}$ zeros out all the rows of a matrix except those corresponding to $A$. Therefore we can split blockwise using the partition $\mathcal{Y}=\bar{A}\cup A$ to get

Note that $K$ can be computed from an eigendecomposition of $L=\sum_{n=1}^{N}\lambda_{n}\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}$ by a simple rescaling of eigenvalues:

Conversely, we can ask when a DPP with marginal kernel $K$ is also an L-ensemble. By inverting Equation (21), we have

and again the computation can be performed by eigendecomposition. However, while the inverse in Equation (21) always exists due to the positive coefficient on the identity matrix, the inverse in Equation (33) may not. In particular, when any eigenvalue of $K$ achieves the upper bound of 1, the DPP is not an L-ensemble. We will see later that the existence of the inverse in Equation (33) is equivalent to $\mathcal{P}$ giving nonzero probability to the empty set. (This is somewhat analogous to the positive probability assumption in the Hammersley-Clifford theorem for Markov random fields.) This is not a major restriction, for two reasons. First, when modeling real data we must typically allocate some nonzero probability for rare or noisy events, so when cardinality is one of the aspects we wish to model, the condition is not unreasonable. Second, we will show in Section 5 how to control the cardinality of samples drawn from the DPP, thus sidestepping the representational limits of L-ensembles.

Modulo the restriction described above, $K$ and $L$ offer alternative representations of DPPs. Under both representations, subsets that have higher diversity, as measured by the corresponding kernel, have higher likelihood. However, while $K$ gives marginal probabilities, L-ensembles directly model the atomic probabilities of observing each subset of $\mathcal{Y}$, which offers an appealing target for optimization. Furthermore, $L$ need only be positive semidefinite, while the eigenvalues of $K$ are bounded above. For these reasons we will focus our modeling efforts on DPPs represented as L-ensembles.

Determinants have an intuitive geometric interpretation. Let $B$ be a $D$ x $N$ matrix such that $L=B^{\top}B$. (Such a $B$ can always be found for $D\leq N$ when $L$ is positive semidefinite.) Denote the columns of $B$ by $B_{i}$ for $i=1,2,\dots,N$. Then:

where the right hand side is the squared $|Y|$-dimensional volume of the parallelepiped spanned by the columns of $B$ corresponding to elements in $Y$.

Intuitively, we can think of the columns of $B$ as feature vectors describing the elements of $\mathcal{Y}$. Then the kernel $L$ measures similarity using dot products between feature vectors, and Equation (34) says that the probability assigned by a DPP to a set $Y$ is related to the volume spanned by its associated feature vectors. This is illustrated in Figure 6.

From this intuition we can verify several important DPP properties. Diverse sets are more probable because their feature vectors are more orthogonal, and hence span larger volumes. Items with parallel feature vectors are selected together with probability zero, since their feature vectors define a degenerate parallelepiped. All else being equal, items with large-magnitude feature vectors are more likely to appear, because they multiply the spanned volumes for sets containing them.

We will revisit these intuitions in Section 3.1, where we decompose the kernel $L$ so as to separately model the direction and magnitude of the vectors $B_{i}$.

3 Properties

In this section we review several useful properties of DPPs.

If $\boldsymbol{Y}$ is distributed as a DPP with marginal kernel $K$, then $\boldsymbol{Y}\cap A$, where $A\subseteq\mathcal{Y}$, is also distributed as a DPP, with marginal kernel $K_{A}$.

If $\boldsymbol{Y}$ is distributed as a DPP with marginal kernel $K$, then $\mathcal{Y}-\boldsymbol{Y}$ is also distributed as a DPP, with marginal kernel $\bar{K}=I-K$. In particular, we have

where $I$ indicates the identity matrix of appropriate size. It may seem counterintuitive that the complement of a diversifying process should also encourage diversity. However, it is easy to see that

where the inequality follows from Equation (7).

If $K\preceq K^{\prime}$, that is, $K^{\prime}-K$ is positive semidefinite, then for all $A\subseteq\mathcal{Y}$ we have

In other words, the DPP defined by $K^{\prime}$ is larger than the one defined by $K$ in the sense that it assigns higher marginal probabilities to every set $A$. An analogous result fails to hold for $L$ due to the normalization constant.

If $K=\gamma K^{\prime}$ for some $0\leq\gamma<1$, then for all $A\subseteq\mathcal{Y}$ we have

It is easy to see that $K$ defines the distribution obtained by taking a random set distributed according to the DPP with marginal kernel $K^{\prime}$, and then independently deleting each of its elements with probability $1-\gamma$.

Note that, by Equation (21), $\frac{\lambda_{1}}{\lambda_{1}+1},\frac{\lambda_{2}}{\lambda_{2}+1},\dots,\frac{\lambda_{N}}{\lambda_{N}+1}$ are the eigenvalues of $K$. Figure 7 shows a plot of the function $f(\lambda)=\frac{\lambda}{\lambda+1}$. It is easy to see from this why the class of L-ensembles does not include DPPs where the empty set has probability zero—at least one of the Bernoulli trials would need to always succeed, and in turn one or more of the eigenvalues of $L$ would be infinite.

In some instances, the sum of Bernoullis may be an appropriate model for uncertain cardinality in real-world data, for instance when identifying objects in images where the number of objects is unknown in advance. In other situations, it may be more practical to fix the cardinality of $\boldsymbol{Y}$ up front, for instance when a set of exactly ten search results is desired, or to replace the sum of Bernoullis with an alternative cardinality model. We show how these goals can be can be achieved in Section 5.

4 Inference

One of the primary advantages of DPPs is that, although the number of possible realizations of $\boldsymbol{Y}$ is exponential in $N$, many types of inference can be performed in polynomial time. In this section we review the inference questions that can (and cannot) be answered efficiently. We also discuss the empirical practicality of the associated computations and algorithms, estimating the largest values of $N$ that can be handled at interactive speeds (within 2–3 seconds) as well as under more generous time and memory constraints. The reference machine used for estimating real-world performance has eight Intel Xeon E5450 3Ghz cores and 32GB of memory.

As we have already seen, the partition function, despite being a sum over $2^{N}$ terms, can be written in closed form as $\det(L+I)$. Determinants of $N\times N$ matrices can be computed through matrix decomposition in $O(N^{3})$ time, or reduced to matrix multiplication for better asymptotic performance. The Coppersmith-Winograd algorithm, for example, can be used to compute determinants in about $O(N^{2.376})$ time. Going forward, we will use $\omega$ to denote the exponent of whatever matrix multiplication algorithm is used.

Practically speaking, modern computers can calculate determinants up to $N\approx 5{,}000$ at interactive speeds, or up to $N\approx 40{,}000$ in about five minutes. When $N$ grows much larger, the memory required simply to store the matrix becomes limiting. (Sparse storage of larger matrices is possible, but computing determinants remains prohibitively expensive unless the level of sparsity is extreme.)

4.2 Marginalization

The marginal probability of any set of items $A$ can be computed using the marginal kernel as in Equation (1). From Equation (35) we can also compute the marginal probability that none of the elements in $A$ appear. (We will see below how marginal probabilities of arbitrary configurations can be computed using conditional DPPs.)

If the DPP is specified as an L-ensemble, then the computational bottleneck for marginalization is the computation of $K$. The dominant operation is the matrix inversion, which requires at least $O(N^{\omega})$ time by reduction to multiplication, or $O(N^{3})$ using Gauss-Jordan elimination or various matrix decompositions, such as the eigendecomposition method proposed in Section 2.2. Since an eigendecomposition of the kernel will be central to sampling, the latter approach is often the most practical when working with DPPs.

Matrices up to $N\approx 2{,}000$ can be inverted at interactive speeds, and problems up to $N\approx 20{,}000$ can be completed in about ten minutes.

4.3 Conditioning

The distribution obtained by conditioning a DPP on the event that none of the elements in a set $A$ appear is easy to compute. For $B\subseteq{\mathcal{Y}}$ not intersecting with $A$ we have

where $L_{\bar{A}}$ is the restriction of $L$ to the rows and columns indexed by elements in $\mathcal{Y}-A$. In other words, the conditional distribution (over subsets of $\mathcal{Y}-A$) is itself a DPP, and its kernel $L_{\bar{A}}$ is obtained by simply dropping the rows and columns of $L$ that correspond to elements in $A$.

We can also condition a DPP on the event that all of the elements in a set $A$ are observed. For $B$ not intersecting with $A$ we have

where $I_{\bar{A}}$ is the matrix with ones in the diagonal entries indexed by elements of $\mathcal{Y}-A$ and zeros everywhere else. Though it is not immediately obvious, Borodin and Rains (2005) showed that this conditional distribution (over subsets of $\mathcal{Y}-A$) is again a DPP, with a kernel given by

(Following the $N\times N$ inversion, the matrix is restricted to rows and columns indexed by elements in $\mathcal{Y}-A$, then inverted again.) It is easy to show that the inverses exist if and only if the probability of $A$ appearing is nonzero.

Combining Equation (46) and Equation (49), we can write the conditional DPP given an arbitrary combination of appearing and non-appearing elements:

Thus, the class of DPPs is closed under most natural conditioning operations.

These formulas also allow us to compute arbitrary marginals. For example, applying Equation (21) to Equation (50) yields the marginal kernel for the conditional DPP given the appearance of $A$:

(Note that $K^{A}$ is indexed by elements of $\mathcal{Y}-A$, so this is only defined when $A$ and $B$ are disjoint.) Using Equation (35) for the complement of a DPP, we can now compute the marginal probability of any partial assignment, i.e.,

Computing conditional DPP kernels in general is asymptotically as expensive as the dominant matrix inversion, although in some cases (conditioning only on non-appearance), the inversion is not necessary. In any case, conditioning is at most a small constant factor more expensive than marginalization.

4.4 Sampling

Algorithm 1, due to Hough et al. (2006), gives an efficient algorithm for sampling a configuration $Y$ from a DPP. The input to the algorithm is an eigendecomposition of the DPP kernel $L$. Note that $\boldsymbol{e}_{i}$ is the $i$th standard basis $N$-vector, which is all zeros except for a one in the $i$th position. We will prove the following theorem.

Let $L=\sum_{n=1}^{N}\lambda_{n}\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}$ be an orthonormal eigendecomposition of a positive semidefinite matrix $L$. Then Algorithm 1 samples $\boldsymbol{Y}\sim\mathcal{P}_{L}$.

Algorithm 1 has two main loops, corresponding to two phases of sampling. In the first phase, a subset of the eigenvectors is selected at random, where the probability of selecting each eigenvector depends on its associated eigenvalue. In the second phase, a sample $Y$ is produced based on the selected vectors. Note that on each iteration of the second loop, the cardinality of $Y$ increases by one and the dimension of $V$ is reduced by one. Since the initial dimension of $V$ is simply the number of selected eigenvectors ($|J|$), Theorem 2.3 has the previously stated corollary that the cardinality of a random sample is distributed as a sum of Bernoulli variables.

To prove Theorem 2.3 we will first show that a DPP can be expressed as a mixture of simpler, elementary DPPs. We will then show that the first phase chooses an elementary DPP according to its mixing coefficient, while the second phase samples from the elementary DPP chosen in phase one.

A DPP is called elementary if every eigenvalue of its marginal kernel is in $\{0,1\}$. We write $\mathcal{P}^{V}$, where $V$ is a set of orthonormal vectors, to denote an elementary DPP with marginal kernel $K^{V}=\sum_{\boldsymbol{v}\in V}\boldsymbol{v}\boldsymbol{v}^{\top}$.

We introduce the term “elementary” here; Hough et al. (2006) refer to elementary DPPs as determinantal projection processes, since $K^{V}$ is an orthonormal projection matrix to the subspace spanned by $V$. Note that, due to Equation (33), elementary DPPs are not generally L-ensembles. We start with a technical lemma.

Let $W_{n}$ for $n=1,2,\dots,N$ be an arbitrary sequence of $k\times k$ rank-one matrices, and let $(W_{n})_{i}$ denote the $ith$ column of $W_{n}$. Let $W_{J}=\sum_{n\in J}W_{n}$. Then

Expanding on the first column of $W_{J}$ using the multilinearity of the determinant,

and, applying the same operation inductively to all columns,

Since $W_{n}$ has rank one, the determinant of any matrix containing two or more columns of $W_{n}$ is zero; thus the terms in the sum vanish unless $n_{1},n_{2},\dots,n_{k}$ are distinct. ∎

A DPP with kernel $L=\sum_{n=1}^{N}\lambda_{n}\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}$ is a mixture of elementary DPPs:

where $V_{J}$ denotes the set $\{\boldsymbol{v}_{n}\}_{n\in J}$.

Consider an arbitrary set $A$, with $k=|A|$. Let $W_{n}=[\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}]_{A}$ for $n=1,2,\dots,N$; note that all of the $W_{n}$ have rank one. From the definition of $K^{V_{J}}$, the mixture distribution on the right hand side of Equation (60) gives the following expression for the marginal probability of $A$:

using the fact that $\det(L+I)=\prod_{n=1}^{N}(\lambda_{n}+1)$. Applying Lemma 57 in reverse and then the definition of $K$ in terms of the eigendecomposition of $L$, we have that the marginal probability of $A$ given by the mixture is

Since the two distributions agree on all marginals, they are equal. ∎

Next, we show that elementary DPPs have fixed cardinality.

If $\boldsymbol{Y}$ is drawn according to an elementary DPP $\mathcal{P}^{V}$, then $|\boldsymbol{Y}|=|V|$ with probability one.

Since $K^{V}$ has rank $|V|$, $\mathcal{P}^{V}(Y\subseteq\boldsymbol{Y})=0$ whenever $|Y|>|V|$, so $|\boldsymbol{Y}|\leq|V|$. But we also have

Thus $|\boldsymbol{Y}|=|V|$ almost surely. ∎

Lemma 2.6 says that the mixture weight of $\mathcal{P}^{V_{J}}$ is given by the product of the eigenvalues $\lambda_{n}$ corresponding to the eigenvectors $\boldsymbol{v}_{n}\in V_{J}$, normalized by $\det(L+I)=\prod_{n=1}^{N}(\lambda_{n}+1)$. This shows that the first loop of Algorithm 1 selects an elementary DPP $\mathcal{P}^{V}$ with probability equal to its mixture component. All that remains is to show that the second loop samples $\boldsymbol{Y}\sim\mathcal{P}^{V}$.

Let $B$ represent the matrix whose rows are the eigenvectors in $V$, so that $K^{V}=B^{\top}B$. Using the geometric interpretation of determinants introduced in Section 2.2.1, $\det(K^{V}_{Y})$ is equal to the squared volume of the parallelepiped spanned by $\{B_{i}\}_{i\in Y}$. Note that since $V$ is an orthonormal set, $B_{i}$ is just the projection of $\boldsymbol{e}_{i}$ onto the subspace spanned by $V$.

Let $k=|V|$. By Lemma 2.7 and symmetry, we can consider without loss of generality a single $Y=\{1,2,\dots,k\}$. Using the fact that any vector both in the span of $V$ and perpendicular to $\boldsymbol{e}_{i}$ is also perpendicular to the projection of $\boldsymbol{e}_{i}$ onto the span of $V$, by the base $\times$ height formula for the volume of a parallelepiped we have

Assume that, as iteration $j$ of the second loop in Algorithm 1 begins, we have already selected $y_{1}=1,y_{2}=2,\dots,y_{j-1}=j-1$. Then $V$ in the algorithm has been updated to an orthonormal basis for the subspace of the original $V$ perpendicular to $\boldsymbol{e}_{1},\dots,\boldsymbol{e}_{j-1}$, and the probability of choosing $y_{j}=j$ is exactly

Therefore, the probability of selecting the sequence $1,2,\dots,k$ is

Since volume is symmetric, the argument holds identically for all of the $k!$ orderings of $Y$. Thus the total probability that Algorithm 1 selects $Y$ is $\det(K^{V}_{Y})$. ∎

Algorithm 1 generates $Y$ in uniformly random order.

To get a feel for the sampling algorithm, it is useful to visualize the distributions used to select $i$ at each iteration, and to see how they are influenced by previously chosen items. Figure 8a shows this progression for a simple DPP where $\mathcal{Y}$ is a finely sampled grid of points in $$, and the kernel is such that points are more similar the closer together they are. Initially, the eigenvectors$V$give rise to a fairly uniform distribution over points in$\mathcal{Y}$, but as each successive point is selected and$V$ is updated, the distribution shifts to avoid points near those already chosen. Figure 8b shows a similar progression for a DPP over points in the unit square.

The sampling algorithm also offers an interesting analogy to clustering. If we think of the eigenvectors of $L$ as representing soft clusters, and the eigenvalues as representing their strengths—the way we do for the eigenvectors and eigenvalues of the Laplacian matrix in spectral clustering—then a DPP can be seen as performing a clustering of the elements in $\mathcal{Y}$, selecting a random subset of clusters based on their strength, and then choosing one element per selected cluster. Of course, the elements are not chosen independently and cannot be identified with specific clusters; instead, the second loop of Algorithm 1 coordinates the choices in a particular way, accounting for overlap between the eigenvectors.

Algorithm 1 runs in time $O(Nk^{3})$, where $k=|V|$ is the number of eigenvectors selected in phase one. The most expensive operation is the $O(Nk^{2})$ Gram-Schmidt orthonormalization required to compute $V_{\perp}$. If $k$ is large, this can be reasonably expensive, but for most applications we do not want high-cardinality DPPs. (And if we want very high-cardinality DPPs, we can potentially save time by using Equation (35) to sample the complement instead.) In practice, the initial eigendecomposition of $L$ is often the computational bottleneck, requiring $O(N^{3})$ time. Modern multi-core machines can compute eigendecompositions up to $N\approx 1{,}000$ at interactive speeds of a few seconds, or larger problems up to $N\approx 10{,}000$ in around ten minutes. In some instances, it may be cheaper to compute only the top $k$ eigenvectors; since phase one tends to choose eigenvectors with large eigenvalues anyway, this can be a reasonable approximation when the kernel is expected to be low rank. Note that when multiple samples are desired, the eigendecomposition needs to be performed only once.

Deshpande and Rademacher (2010) recently proposed a $(1-\epsilon)$-approximate algorithm for sampling that runs in time $O(N^{2}\log N\frac{k^{2}}{\epsilon^{2}}+N\log^{\omega}N\frac{k^{2\omega+1}}{\epsilon^{2\omega}}\log(\frac{k}{\epsilon}\log N))$ when $L$ is already decomposed as a Gram matrix, $L=B^{\top}B$. When $B$ is known but an eigendecomposition is not (and $N$ is sufficiently large), this may be significantly faster than the exact algorithm.

4.5 Finding the mode

Finding the mode of a DPP—that is, finding the set $Y\subseteq\mathcal{Y}$ that maximizes $\mathcal{P}_{L}(Y)$—is NP-hard. In conditional models, this problem is sometimes referred to as maximum a posteriori (or MAP) inference, and it is also NP-hard for most general structured models such as Markov random fields. Hardness was first shown for DPPs by Ko et al. (1995), who studied the closely-related maximum entropy sampling problem: the entropy of a set of jointly Gaussian random variables is given (up to constants) by the log-determinant of their covariance matrix; thus finding the maximum entropy subset of those variables requires finding the principal covariance submatrix with maximum determinant. Here, we adapt the argument of Çivril and Magdon-Ismail (2009), who studied the problem of finding maximum-volume submatrices.

Let dpp-mode be the optimization problem of finding, for a positive semidefinite $N\times N$ input matrix $L$ indexed by elements of $\mathcal{Y}$, the maximum value of $\det(L_{Y})$ over all $Y\subseteq\mathcal{Y}$. dpp-mode is NP-hard, and furthermore it is NP-hard even to approximate dpp-mode to a factor of $\frac{8}{9}+\epsilon$.

We reduce from exact 3-cover (X3C). An instance of X3C is a set $S$ and a collection $C$ of three-element subsets of $S$; the problem is to decide whether there is a sub-collection $C^{\prime}\subseteq C$ such that every element of $S$ appears exactly once in $C^{\prime}$ (that is, $C^{\prime}$ is an exact 3-cover). X3C is known to be NP-complete.

The reduction is as follows. Let $\mathcal{Y}=\{1,2,\dots,|C|\}$, and let $B$ be an $|S|\times|C|$ matrix where $B_{si}=\frac{1}{\sqrt{3}}$ if $C_{i}$ contains $s\in S$ and zero otherwise. Define $L=\gamma B^{\top}B$, where $1<\gamma\leq\frac{9}{8}$. Note that the diagonal of $L$ is constant and equal to $\gamma$, and an off-diagonal entry $L_{ij}$ is zero if and only if $C_{i}$ and $C_{j}$ do not intersect. $L$ is positive semidefinite by construction, and the reduction requires only polynomial time. Let $k=\frac{|S|}{3}$. We will show that the maximum value of $\det(L_{Y})$ is greater than $\gamma^{k-1}$ if and only if $C$ contains an exact 3-cover of $S$.

$(\leftarrow)$ If $C^{\prime}\subseteq C$ is an exact 3-cover of $S$, then it must contain exactly $k$ 3-sets. Letting $Y$ be the set of indices in $C^{\prime}$, we have $L_{Y}=\gamma I$, and thus its determinant is $\gamma^{k}>\gamma^{k-1}$.

$(\rightarrow)$ Suppose there is no 3-cover of $S$ in $C$. Let $Y$ be an arbitrary subset of $\mathcal{Y}$. If $|Y|<k$, then

Now suppose $|Y|\geq k$, and assume without loss of generality that $Y=\{1,2,\dots,|Y|\}$. We have $L_{Y}=\gamma B_{Y}^{\top}B_{Y}$, and

Note that, since the columns of $B$ are normalized, each term in the product is at most one. Furthermore, at least $|Y|-k+1$ of the terms must be strictly less than one, because otherwise there would be $k$ orthogonal columns, which would correspond to a 3-cover. By the construction of $B$, if two columns $B_{i}$ and $B_{j}$ are not orthogonal then $C_{i}$ and $C_{j}$ overlap in at least one of three elements, so we have

We have shown that the existence of a 3-cover implies that the optimal value of dpp-mode is at least $\gamma^{k}$, while the optimal value cannot be more than $\gamma^{k-1}$ if there is no 3-cover. Thus any algorithm that can approximate dpp-mode to better than a factor of $\frac{1}{\gamma}$ can be used to solve X3C in polynomial time. We can choose $\gamma=\frac{9}{8}$ to show that an approximation ratio of $\frac{8}{9}+\epsilon$ is NP-hard. ∎

Since there are only $|C|$ possible cardinalities for $Y$, Theorem 2.9 shows that dpp-mode is NP-hard even under cardinality constraints.

(Ko et al., 1995) propose an exact, exponential branch-and-bound algorithm for finding the mode using greedy heuristics to build candidate sets; they tested their algorithm on problems up to $N=75$, successfully finding optimal solutions in up to about an hour. Modern computers are likely a few orders of magnitude faster; however, this algorithm is still probably impractical for applications with large $N$. Çivril and Magdon-Ismail (2009) propose an efficient greedy algorithm for finding a set of size $k$, and prove that it achieves an approximation ratio of $O(\frac{1}{k!})$. While this guarantee is relatively poor for all but very small $k$, in practice the results may be useful nonetheless.

$\mathcal{P}_{L}$ is log-submodular; that is,

whenever $Y\subseteq Y^{\prime}\subseteq\mathcal{Y}-\{i\}$. Intuitively, adding elements to $Y$ yields diminishing returns as $Y$ gets larger. (This is easy to show by a volume argument.) Submodular functions can be minimized in polynomial time (Schrijver, 2000), and many results exist for approximately maximizing monotone submodular functions, which have the special property that supersets always have higher function values than their subsets (Nemhauser et al., 1978; Fisher et al., 1978; Feige, 1998). In Section 4.2.1 we will discuss how these kinds of greedy algorithms can be adapted for DPPs. However, in general $\mathcal{P}_{L}$ is highly non-monotone, since the addition of even a single element can decrease the probability to zero.

Recently, Feige et al. (2007) showed that even non-monotone submodular functions can be approximately maximized in polynomial time using a local search algorithm, and a growing body of research has focused on extending this result in a variety of ways (Lee et al., 2009; Gharan and Vondrák, 2011; Vondrák et al., 2011; Feldman et al., 2011a, b; Chekuri et al., 2011). In our recent work we showed how the computational structure of DPPs gives rise to a particularly efficient variant of these methods (Kulesza et al., 2012).

5 Related processes

Historically, a wide variety of point process models have been proposed and applied to applications involving diverse subsets, particularly in settings where the items can be seen as points in a physical space and diversity is taken to mean some sort of “spreading” behavior. However, DPPs are essentially unique among this class in having efficient and exact algorithms for probabilistic inference, which is why they are particularly appealing models for machine learning applications. In this section we briefly survey the wider world of point processes and discuss the computational properties of alternative models; we will focus on point processes that lead to what is variously described as diversity, repulsion, (over)dispersion, regularity, order, and inhibition.

Perhaps the most fundamental point process is the Poisson point process, which is depicted on the right side of Figure 4 (Daley and Vere-Jones, 2003). While defined for continuous $\mathcal{Y}$, in the discrete setting the Poisson point process can be simulated by flipping a coin independently for each item, and including those items for which the coin comes up heads. Formally,

where $p_{i}\in$ is the bias of the $i$th coin. The process is called stationary when $p_{i}$ does not depend on $i$; in a spatial setting this means that no region has higher density than any other region.

A random set $\boldsymbol{Y}$ distributed as a Poisson point process has the property that whenever $A$ and $B$ are disjoint subsets of $\mathcal{Y}$, the random variables $\boldsymbol{Y}\cap A$ and $\boldsymbol{Y}\cap B$ are independent; that is, the points in $\boldsymbol{Y}$ are not correlated. It is easy to see that DPPs generalize Poisson point processes by choosing the marginal kernel $K$ with $K_{ii}=p_{i}$ and $K_{ij}=0,i\neq j$. This implies that inference for Poisson point processes is at least as efficient as for DPPs; in fact, it is more efficient, since for instance it is easy to compute the most likely configuration. However, since Poisson point processes do not model correlations between variables, they are rather uninteresting for most real-world applications.

Addressing this weakness, various procedural modifications of the Poisson process have been proposed in order to introduce correlations between items. While such constructions can be simple and intuitive, leading to straightforward sampling algorithms, they tend to make general statistical inference difficult.

Matérn (1960, 1986) proposed a set of techniques for thinning Poisson point processes in order to induce a type of repulsion when the items are embedded in a Euclidean space. The Type I process is obtained from a Poisson set $\boldsymbol{Y}$ by removing all items in $\boldsymbol{Y}$ that lie within some radius of another item in $\boldsymbol{Y}$. That is, if two items are close to each other, they are both removed; as a result all items in the final process are spaced at least a fixed distance apart. The Type II Matérn repulsive process, designed to achieve the same minimum distance property while keeping more items, begins by independently assigning each item in $\boldsymbol{Y}$ a uniformly random “time” in $$. Then, any item within a given radius of a point having a smaller time value is removed. Under this construction, when two items are close to each other only the later one is removed. Still, an item may be removed due to its proximity with an earlier item that was itself removed. This leads to the Type III process, which proceeds dynamically, eliminating items in time order whenever an earlier point which has not been removed lies within the radius.

Inference for the Matérn processes is computationally daunting. First and second order moments can be computed for Types I and II, but in those cases computing the likelihood of a set $Y$ is seemingly intractable (Møller et al., 2010). Recent work by Huber and Wolpert (2009) shows that it is possible to compute likelihood for certain restricted Type III processes, but computing moments cannot be done in closed form. In the general case, likelihood for Type III processes must be estimated using an expensive Markov chain Monte Carlo algorithm.

The Matérn processes are called “hard-core” because they strictly enforce a minimum radius between selected items. While this property leads to one kind of diversity, it is somewhat limited, and due to the procedural definition it is difficult to characterize the side effects of the thinning process in a general way. Stoyan and Stoyan (1985) considered an extension where the radius is itself chosen randomly, which may be more natural for certain settings, but it does not alleviate the computational issues.

The Matérn repulsive processes are related in spirit to the random sequential adsorption (RSA) model, which has been used in physics and chemistry to model particles that bind to two-dimensional surfaces, e.g., proteins on a cell membrane (Tanemura, 1979; Finegold and Donnell, 1979; Feder, 1980; Swendsen, 1981; Hinrichsen et al., 1986; Ramsden, 1993). RSA is described generatively as follows. Initially, the surface is empty; iteratively, particles arrive and bind uniformly at random to a location from among all locations that are not within a given radius of any previously bound particle. When no such locations remain (the “jamming limit”), the process is complete.

Like the Matérn processes, RSA is a hard-core model, designed primarily to capture packing distributions, with much of the theoretical analysis focused on the achievable density. If the set of locations is further restricted at each step to those found in an initially selected Poisson set $\boldsymbol{Y}$, then it is equivalent to a Matérn Type III process (Huber and Wolpert, 2009); it therefore shares the same computational burdens.

5.2 Gibbs and Markov point processes

While manipulating the Poisson process procedurally has some intuitive appeal, it seems plausible that a more holistically-defined process might be easier to work with, both analytically and algorithmically. The Gibbs point process provides such an approach, offering a general framework for incorporating correlations among selected items (Preston, 1976; Ripley and Kelly, 1977; Ripley, 1991; Van Lieshout, 2000; Møller and Waagepetersen, 2004, 2007; Daley and Vere-Jones, 2008). The Gibbs probability of a set $Y$ is given by

where $U$ is an energy function. Of course, this definition is fully general without further constraints on $U$. A typical assumption is that $U$ decomposes over subsets of items in $Y$; for instance

for some small constant order $k$ and potential functions $\psi$. In practice, the most common case is $k=2$, which is sometimes called a pairwise interaction point process (Diggle et al., 1987):

In spatial settings, a Gibbs point process whose potential functions are identically 1 whenever their input arguments do not lie within a ball of fixed radius—that is, whose energy function can be decomposed into only local terms—is called a Markov point process. A number of specific Markov point processes have become well-known.

Strauss (1975) introduced a simple pairwise Markov point process for spatial data in which the potential function $\psi_{2}(i,j)$ is piecewise constant, taking the value 1 whenever $i$ and $j$ are at least a fixed radius apart, and the constant value $\gamma$ otherwise. When $\gamma>1$, the resulting process prefers clustered items. (Note that $\gamma>1$ is only possible in the discrete case; in the continuous setting the distribution becomes non-integrable.) We are more interested in the case $0<\gamma<1$, where configurations in which selected items are near one another are discounted. When $\gamma=0$, the resulting process becomes hard-core, but in general the Strauss process is “soft-core”, preferring but not requiring diversity.

The Strauss process is typical of pairwise Markov processes in that its potential function $\psi_{2}(i,j)=\psi(|i-j|)$ depends only on the distance between its arguments. A variety of alternative definitions for $\psi(\cdot)$ have been proposed (Ripley and Kelly, 1977; Ogata and Tanemura, 1984). For instance,

where $\sigma$ controls the degree of repulsion in each case. Each definition leads to a point process with a slightly different concept of diversity.

Baddeley and Van Lieshout (1995) proposed a non-pairwise spatial Markov point process called the area-interaction model, where $U(Y)$ is given by $\log\gamma$ times the total area contained in the union of discs of fixed radius centered at all of the items in $Y$. When $\gamma>1$, we have $\log\gamma>0$ and the process prefers sets whose discs cover as little area as possible, i.e., whose items are clustered. When $0<\gamma<1$, $\log\gamma$ becomes negative, so the process prefers “diverse” sets covering as much area as possible.

If none of the selected items fall within twice the disc radius of each other, then $\exp(-U(Y))$ can be decomposed into potential functions over single items, since the total area is simply the sum of the individual discs. Similarly, if each disc intersects with at most one other disc, the area-interaction process can be written as a pairwise interaction model. However, in general, an unbounded number of items might appear in a given disc; as a result the area-interaction process is an infinite-order Gibbs process. Since items only interact when they are near one another, however, local potential functions are sufficient and the process is Markov.

Markov point processes have many intuitive properties. In fact, it is not difficult to see that, for discrete ground sets $\mathcal{Y}$, the Markov point process is equivalent to a Markov random field (MRF) on binary variables corresponding to the elements of $\mathcal{Y}$. In Section 3.2.2 we will return to this equivalence in order to discuss the relative expressive possibilities of DPPs and MRFs. For now, however, we simply observe that, as for MRFs with negative correlations, repulsive Markov point processes are computationally intractable. Even computing the normalizing constant for Equation (84) is NP-hard in the cases outlined above (Daley and Vere-Jones, 2003; Møller and Waagepetersen, 2004).

On the other hand, quite a bit of attention has been paid to approximate inference algorithms for Markov point processes, employing pseudolikelihood (Besag, 1977; Besag et al., 1982; Jensen and Moller, 1991; Ripley, 1991), Markov chain Monte Carlo methods (Ripley and Kelly, 1977; Besag and Green, 1993; Häggström et al., 1999; Berthelsen and Møller, 2006), and other approximations (Ogata and Tanemura, 1985; Diggle et al., 1994). Nonetheless, in general these methods are slow and/or inexact, and closed-form expressions for moments and densities rarely exist (Møller and Waagepetersen, 2007). In this sense the DPP is unique.

5.3 Generalizations of determinants

The determinant of a $k\times k$ matrix $K$ can be written as a polynomial of degree $k$ in the entries of $K$; in particular,

Computationally, obtaining the permanent of a matrix is #P-complete (Valiant, 1979), making the permanental process difficult to work with in practice. Complexity results for immanants are less definitive, with certain classes of immanants apparently hard to compute (Bürgisser, 2000; Brylinski and Brylinski, 2003), while some upper bounds on complexity are known (Hartmann, 1985; Barvinok, 1990), and at least one non-trivial case is efficiently computable (Grone and Merris, 1984). It is not clear whether the latter result provides enough leverage to perform inference beyond computing marginals.

A third possible generalization of Equation (90) is the hyperdeterminant originally proposed by Cayley (1843) and discussed in the context of point processes by Evans and Gottlieb (2009). Whereas the standard determinant operates on a two-dimensional matrix with entries indexed by pairs of items, the hyperdeterminant operates on higher-dimensional kernel matrices indexed by sets of items. The hyperdeterminant potentially offers additional modeling power, and Evans and Gottlieb (2009) show that some useful properties of DPPs are preserved in this setting. However, so far relatively little is known about these processes.

5.4 Quasirandom processes

Monte Carlo methods rely on draws of random points in order to approximate quantities of interest; randomness guarantees that, regardless of the function being studied, the estimates will be accurate in expectation and converge in the limit. However, in practice we get to observe only a finite set of values drawn from the random source. If, by chance, this set is “bad”, the resulting estimate may be poor. This concern has led to the development of so-called quasirandom sets, which are in fact deterministically generated, but can be substituted for random sets in some instances to obtain improved convergence guarantees (Niederreiter, 1992; Sobol, 1998).

In contrast with pseudorandom generators, which attempt to mimic randomness by satisfying statistical tests that ensure unpredictability, quasirandom sets are not designed to appear random, and their elements are not (even approximately) independent. Instead, they are designed to have low discrepancy; roughly speaking, low-discrepancy sets are “diverse” in that they cover the sample space evenly. Consider a finite subset $Y$ of $^{D}$, with elements $\boldsymbol{x}^{(i)}=(x^{(i)}_{1},x^{(i)}_{2},\dots,x^{(i)}_{D})$ for $i=1,2,\dots,k$. Let $S_{\boldsymbol{x}}=[0,x_{1})\times[0,x_{2})\times\cdots\times[0,x_{D})$ denote the box defined by the origin and the point $\boldsymbol{x}$. The discrepancy of $Y$ is defined as follows.

This deterministic uniformity property makes quasirandom sets useful for Monte Carlo estimation via (among other results) the Koksma-Hlawka inequality (Hlawka, 1961; Niederreiter, 1992). For a function $f$ with bounded variation $V(f)$ on the unit cube, the inequality states that

Thus, low-discrepancy sets lead to accurate quasi-Monte Carlo estimates. In contrast to typical Monte Carlo guarantees, the Koksma-Hlawka inequality is deterministic. Moreover, since the rate of convergence for standard stochastic Monte Carlo methods is $k^{-1/2}$, this result is an (asymptotic) improvement when the discrepancy diminishes faster than $k^{-1/2}$.

In fact, it is possible to construct quasirandom sequences where the discrepancy of the first $k$ elements is $O((\log k)^{D}/k)$; the first such sequence was proposed by Halton (1960). The Sobol sequence (Sobol, 1967), introduced later, offers improved uniformity properties and can be generated efficiently (Bratley and Fox, 1988).

It seems plausible that, due to their uniformity characteristics, low-discrepancy sets could be used as computationally efficient but non-probabilistic tools for working with data exhibiting diversity. An algorithm generating quasirandom sets could be seen as an efficient prediction procedure if made to depend somehow on input data and a set of learned parameters. However, to our knowledge no work has yet addressed this possibility.

Representation and algorithms

Determinantal point processes come with a deep and beautiful theory, and, as we have seen, exactly characterize many theoretical processes. However, they are also promising models for real-world data that exhibit diversity, and we are interested in making such applications as intuitive, practical, and computationally efficient as possible. In this section, we present a variety of fundamental techniques and algorithms that serve these goals and form the basis of the extensions we discuss later.

We begin by describing a decomposition of the DPP kernel that offers an intuitive tradeoff between a unary model of quality over the items in the ground set and a global model of diversity. The geometric intuitions from Section 2 extend naturally to this decomposition. Splitting the model into quality and diversity components then allows us to make a comparative study of expressiveness—that is, the range of distributions that the model can describe. We compare the expressive powers of DPPs and negative-interaction Markov random fields, showing that the two models are incomparable in general but exhibit qualitatively similar characteristics, despite the computational advantages offered by DPPs.

Next, we turn to the challenges imposed by large data sets, which are common in practice. We first address the case where $N$, the number of items in the ground set, is very large. In this setting, the super-linear number of operations required for most DPP inference algorithms can be prohibitively expensive. However, by introducing a dual representation of a DPP we show that efficient DPP inference remains possible when the kernel is low-rank. When the kernel is not low-rank, we prove that a simple approximation based on random projections dramatically speeds inference while guaranteeing that the deviation from the original distribution is bounded. These techniques will be especially useful in Section 6, when we consider exponentially large $N$.

Finally, we discuss some alternative formulas for the likelihood of a set $Y$ in terms of the marginal kernel $K$. Compared to the L-ensemble formula in Equation (19), these may be analytically more convenient, since they do not involve ratios or arbitrary principal minors.

An important practical concern for modeling is interpretability; that is, practitioners should be able to understand the parameters of the model in an intuitive way. While the entries of the DPP kernel are not totally opaque in that they can be seen as measures of similarity—reflecting our primary qualitative characterization of DPPs as diversifying processes—in most practical situations we want diversity to be balanced against some underlying preferences for different items in $\mathcal{Y}$. In this section, we propose a decomposition of the DPP that more directly illustrates the tension between diversity and a per-item measure of quality.

This decomposition of $L$ has two main advantages. First, it implicitly enforces the constraint that $L$ must be positive semidefinite, which can potentially simplify learning (see Section 4). Second, it allows us to independently model quality and diversity, and then combine them into a unified model. In particular, we have:

where the first term increases with the quality of the selected items, and the second term increases with the diversity of the selected items. We will refer to $q$ as the quality model and $S$ or $\phi$ as the diversity model. Without the diversity model, we would choose high-quality items, but we would tend to choose similar high-quality items over and over. Without the quality model, we would get a very diverse set, but we might fail to include the most important items in $\mathcal{Y}$, focusing instead on low-quality outliers. By combining the two models we can achieve a more balanced result.

Returning to the geometric intuitions from Section 2.2.1, the determinant of $L_{Y}$ is equal to the squared volume of the parallelepiped spanned by the vectors $q_{i}\phi_{i}$ for $i\in Y$. The magnitude of the vector representing item $i$ is $q_{i}$, and its direction is $\phi_{i}$. Figure 9 (reproduced from the previous section) now makes clear how DPPs decomposed in this way naturally balance the two objectives of high quality and high diversity. Going forward, we will nearly always assume that our models are decomposed into quality and diversity components. This provides us not only with a natural and intuitive setup for real-world applications, but also a useful perspective for comparing DPPs with existing models, which we turn to next.

2 Expressive power

Many probabilistic models are known and widely used within the machine learning community. A natural question, therefore, is what advantages DPPs offer that standard models do not. We have seen already how a large variety of inference tasks, like sampling and conditioning, can be performed efficiently for DPPs; however efficiency is essentially a prerequisite for any practical model. What makes DPPs particularly unique is the marriage of these computational advantages with the ability to express global, negative interactions between modeling variables; this repulsive domain is notoriously intractable using traditional approaches like graphical models (Murphy et al., 1999; Boros and Hammer, 2002; Ishikawa, 2003; Taskar et al., 2004; Yanover and Weiss, 2002; Yanover et al., 2006; Kulesza and Pereira, 2008). In this section we elaborate on the expressive powers of DPPs and compare them with those of Markov random fields, which we take as representative graphical models.

A Markov random field (MRF) is an undirected graphical model defined by a graph $G$ whose nodes $1,2,\dots,N$ represent random variables. For our purposes, we will consider binary MRFs, where each output variable takes a value from $\{0,1\}$. We use $y_{i}$ to denote a value of the $i$th output variable, bold $\boldsymbol{y}_{c}$ to denote an assignment to a set of variables $c$, and $\boldsymbol{y}$ for an assignment to all of the output variables. The graph edges $E$ encode direct dependence relationships between variables; for example, there might be edges between similar elements $i$ and $j$ to represent the fact that they tend not to co-occur. MRFs are often referred to as conditional random fields when they are parameterized to depend on input data, and especially when $G$ is a chain (Lafferty et al., 2001).

An MRF defines a joint probability distribution over the output variables that factorizes across the cliques $\mathcal{C}$ of $G$:

Here each $\psi_{c}$ is a potential function that assigns a nonnegative value to every possible assignment $\boldsymbol{y}_{c}$ of the clique $c$, and $Z$ is the normalization constant $\sum_{\boldsymbol{y}^{\prime}}\prod_{c\in\mathcal{C}}\psi_{c}(\boldsymbol{y}^{\prime}_{c})$. Note that, for a binary MRF, we can think of $\boldsymbol{y}$ as the characteristic vector for a subset $Y$ of $\mathcal{Y}=\{1,2,\dots,N\}$. Then the MRF is equivalently the distribution of a random subset $\boldsymbol{Y}$, where $\mathcal{P}(\boldsymbol{Y}=Y)$ is equivalent to $\mathcal{P}(\boldsymbol{y})$.

The Hammersley-Clifford theorem states that $\mathcal{P}(\boldsymbol{y})$ defined in Equation (96) is always Markov with respect to $G$; that is, each variable is conditionally independent of all other variables given its neighbors in $G$. The converse also holds: any distribution that is Markov with respect to $G$, as long as it is strictly positive, can be decomposed over the cliques of $G$ as in Equation (96) (Grimmett, 1973). MRFs therefore offer an intuitive way to model problem structure. Given domain knowledge about the nature of the ways in which outputs interact, a practitioner can construct a graph that encodes a precise set of conditional independence relations. (Because the number of unique assignments to a clique $c$ is exponential in $|c|$, computational constraints generally limit us to small cliques.)

For comparison with DPPs, we will focus on pairwise MRFs, where the largest cliques with interesting potential functions are the edges; that is, $\psi_{c}(\boldsymbol{y}_{c})=1$ for all cliques $c$ where $|c|>2$. The pairwise distribution is

We refer to $\psi_{i}(y_{i})$ as node potentials, and $\psi_{ij}(y_{i},y_{j})$ as edge potentials.

MRFs are very general models—in fact, if the cliques are unbounded in size, they are fully general—but inference is only tractable in certain special cases. Cooper (1990) showed that general probabilistic inference (conditioning and marginalization) in MRFs is NP-hard, and this was later extended by Dagum and Luby (1993), who showed that inference is NP-hard even to approximate. Shimony (1994) proved that the MAP inference problem (finding the mode of an MRF) is also NP-hard, and Abdelbar and Hedetniemi (1998) showed that the MAP problem is likewise hard to approximate. In contrast, we showed in Section 2 that DPPs offer efficient exact probabilistic inference; furthermore, although the MAP problem for DPPs is NP-hard, it can be approximated to a constant factor under cardinality constraints in polynomial time.

The first tractable subclass of MRFs was identified by Pearl (1982), who showed that belief propagation can be used to perform inference in polynomial time when $G$ is a tree. More recently, certain types of inference in binary MRFs with associative (or submodular) potentials $\psi$ have been shown to be tractable (Boros and Hammer, 2002; Taskar et al., 2004; Kolmogorov and Zabih, 2004). Inference in non-binary associative MRFs is NP-hard, but can be efficiently approximated to a constant factor depending on the size of the largest clique (Taskar et al., 2004). Intuitively, an edge potential is called associative if it encourages the endpoint nodes take the same value (e.g., to be both in or both out of the set $Y$). More formally, associative potentials are at least one whenever the variables they depend on are all equal, and exactly one otherwise.

We can rewrite the pairwise, binary MRF of Equation (97) in a canonical log-linear form:

Here we have eliminated redundancies by forcing $\psi_{i}(0)=1$, $\psi_{ij}(0,0)=\psi_{ij}(0,1)=\psi_{ij}(1,0)=1$, and setting $w_{i}=\log\psi_{i}(1)$, $w_{ij}=\log\psi_{ij}(1,1)$. This parameterization is sometimes called the fully visible Boltzmann machine. Under this representation, the MRF is associative whenever $w_{ij}\geq 0$ for all $ij\in E$.

We have seen that inference in MRFs is tractable when we restrict the graph to a tree or require the potentials to encourage agreement. However, the repulsive potentials necessary to build MRFs exhibiting diversity are the opposite of associative potentials (since they imply $w_{ij}<0$), and lead to intractable inference for general graphs. Indeed, such negative potentials can create “frustrated cycles”, which have been used both as illustrations of common MRF inference algorithm failures (Kulesza and Pereira, 2008) and as targets for improving those algorithms (Sontag and Jaakkola, 2007). A wide array of (informally) approximate inference algorithms have been proposed to mitigate tractability problems, but none to our knowledge effectively and reliably handles the case where potentials exhibit strong repulsion.

2.2 Comparing DPPs and MRFs

Despite the computational issues outlined in the previous section, MRFs are popular models and, importantly, intuitive for practitioners, both because they are familiar and because their potential functions directly model simple, local relationships. We argue that DPPs have a similarly intuitive interpretation using the decomposition in Section 3.1. Here, we compare the distributions realizable by DPPs and MRFs to see whether the tractability of DPPs comes at a large expressive cost.

Consider a DPP over $\mathcal{Y}=\{1,2,\dots,N\}$ with $N\times N$ kernel matrix $L$ decomposed as in Section 3.1; we have

Both of these models can capture negative correlations between indicator variables $y_{i}$. Both models also have $\frac{N(N+1)}{2}$ parameters: the DPP has quality scores $q_{i}$ and similarity measures $S_{ij}$, and the MRF has node log-potentials $w_{i}$ and edge log-potentials $w_{ij}$. The key representational difference is that, while $w_{ij}$ are individually constrained to be nonpositive, the positive semidefinite constraint on the DPP kernel is global. One consequence is that, as a side effect, the MRF can actually capture certain limited positive correlations; for example, a 3-node MRF with $w_{12},w_{13}<0$ and $w_{23}=0$ induces a positive correlation between nodes two and three by virtue of their mutual disagreement with node one. As we have seen in Section 2, the semidefinite constraint prevents the DPP from forming any positive correlations.

More generally, semidefiniteness means that the DPP diversity feature vectors must satisfy the triangle inequality, leading to

for all $i,j,k\in\mathcal{Y}$ since $\|\phi_{i}-\phi_{j}\|\propto\sqrt{1-S_{ij}}$. The similarity measure therefore obeys a type of transitivity, with large $S_{ij}$ and $S_{jk}$ implying large $S_{ik}$.

Equation (101) is not, by itself, sufficient to guarantee that $L$ is positive semidefinite, since $S$ must also be realizable using unit length feature vectors. However, rather than trying to develop further intuition algebraically, we turn to visualization. While it is difficult to depict the feasible distributions of DPPs and MRFs in high dimensions, we can get a feel for their differences even with a small number of elements $N$.

When $N=2$, it is easy to show that the two models are equivalent, in the sense that they can both represent any distribution with negative correlations:

When $N=3$, differences start to become apparent. In this setting both models have six parameters: for the DPP they are $(q_{1},q_{2},q_{3},S_{12},S_{13},S_{23})$, and for the MRF they are $(w_{1},w_{2},w_{3},w_{12},w_{13},w_{23})$. To place the two models on equal footing, we represent each as the product of unnormalized per-item potentials $\psi_{1},\psi_{2},\psi_{3}$ and a single unnormalized ternary potential $\psi_{123}$. This representation corresponds to a factor graph with three nodes and a single, ternary factor (see Figure 10). The probability of a set $Y$ is then given by

Figure 11a depicts four such slices of the realizable DPP distributions, and Figure 11b shows the same slices of the realizable MRF distributions. Points closer to the origin (on the lower left) correspond to “more repulsive” distributions, where the three elements of $\mathcal{Y}$ are less likely to appear together. When $\psi_{123}(1,1,1)$ is large (gray surfaces), negative correlations are weak and the two models give rise to qualitatively similar distributions. As the value of the $\psi_{123}(1,1,1)$ shrinks to zero (red surfaces), the two models become quite different. MRFs, for example, can describe distributions where the first item is strongly anti-correlated with both of the others, but the second and third are not anti-correlated with each other. The transitive nature of the DPP makes this impossible.

To improve visibility, we have constrained $S_{12},S_{13},S_{23}\geq 0$ in Figure 11a. Figure 11c shows a single slice without this constraint; allowing negative similarity makes it possible to achieve strong three-way repulsion with less pairwise repulsion, closing the surface away from the origin. The corresponding MRF slice is shown in Figure 11d, and the two are overlaid in Figure 11e and Figure 11f. Even though there are relatively strong interactions in these plots ($\psi_{123}(1,1,1)=0.1$), the models remain roughly comparable in terms of the distributions they can express.

As $N$ gets larger, we conjecture that the story is essentially the same. DPPs are primarily constrained by a notion of transitivity on the similarity measure; thus it would be difficult to use a DPP to model, for example, data where items repel “distant” items rather than similar items—if $i$ is far from $j$ and $j$ is far from $k$ we cannot necessarily conclude that $i$ is far from $k$. One way of looking at this is that repulsion of distant items induces positive correlations between the selected items, which a DPP cannot represent.

MRFs, on the other hand, are constrained by their local nature and do not effectively model data that are “globally” diverse. For instance, a pairwise MRF we cannot exclude a set of three or more items without excluding some pair of those items. More generally, an MRF assumes that repulsion does not depend on (too much) context, so it cannot express that, say, there can be only a certain number of selected items overall. The DPP can naturally implement this kind of restriction though the rank of the kernel.

3 Dual representation

The standard inference algorithms for DPPs rely on manipulating the kernel $L$ through inversion, eigendecomposition, and so on. However, in situations where $N$ is large we may not be able to work efficiently with $L$—in some cases we may not even have the memory to write it down. In this section, instead, we develop a dual representation of a DPP that shares many important properties with the kernel $L$ but is often much smaller. Afterwards, we will show how this dual representation can be applied to perform efficient inference.

Let $B$ be the $D\times N$ matrix whose columns are given by $B_{i}=q_{i}\phi_{i}$, so that $L=B^{\top}B$. Consider now the matrix

By construction, $C$ is symmetric and positive semidefinite. In contrast to $L$, which is too expensive to work with when $N$ is large, $C$ is only $D\times D$, where $D$ is the dimension of the diversity feature function $\phi$. In many practical situations, $D$ is fixed by the model designer, while $N$ may grow without bound as new items become available; for instance, a search engine may continually add to its database of links. Furthermore, we have the following result.

The nonzero eigenvalues of $C$ and $L$ are identical, and the corresponding eigenvectors are related by the matrix $B$. That is,

is an eigendecomposition of $C$ if and only if

In the forward direction, we assume that $\{(\lambda_{n},\boldsymbol{\hat{v}}_{n})\}_{n=1}^{D}$ is an eigendecomposition of $C$. We have

since $\boldsymbol{\hat{v}}_{n}$ are orthonormal by assumption. Furthermore, for any $n$ we have

using the fact that $C\boldsymbol{\hat{v}}_{n}=\lambda_{n}\boldsymbol{\hat{v}}_{n}$ since $\boldsymbol{\hat{v}}_{n}$ is an eigenvector of $C$. Finally, for any distinct $1\leq a,b\leq D$, we have

Thus $\left\{\left(\lambda_{n},\frac{1}{\sqrt{\lambda_{n}}}B^{\top}\boldsymbol{\hat{v}}_{n}\right)\right\}_{n=1}^{D}$ is an eigendecomposition of $L$. In the other direction, an analogous argument applies once we observe that, since $L=B^{\top}B$, $L$ has rank at most $D$ and therefore at most $D$ nonzero eigenvalues. ∎

Proposition 3.1 shows that $C$ contains quite a bit of information about $L$. In fact, $C$ is sufficient to perform nearly all forms of DPP inference efficiently, including normalization and marginalization in constant time with respect to $N$, and sampling in time linear in $N$.

Recall that the normalization constant for a DPP is given by $\det(L+I)$. If $\lambda_{1},\lambda_{2},\dots,\lambda_{N}$ are the eigenvalues of $L$, then the normalization constant is equal to $\prod_{n=1}^{N}(\lambda_{n}+1)$, since the determinant is the product of the eigenvalues of its argument. By Proposition 3.1, the nonzero eigenvalues of $L$ are also the eigenvalues of the dual representation $C$. Thus, we have

Computing the determinant of $C+I$ requires $O(D^{\omega})$ time.

3.2 Marginalization

Standard DPP marginalization makes use of the marginal kernel $K$, which is of course as large as $L$. However, the dual representation $C$ can be used to compute the entries of $K$. We first eigendecompose the dual representation as $C=\sum_{n=1}^{D}\lambda_{n}\boldsymbol{\hat{v}}_{n}\boldsymbol{\hat{v}}_{n}^{\top}$, which requires $O(D^{\omega})$ time. Then, we can use the definition of $K$ in terms of the eigendecomposition of $L$ as well as Proposition 3.1 to compute

That is, the diagonal entries of $K$ are computable from the dot products between the diversity features $\phi_{i}$ and the eigenvectors of $C$; we can therefore compute the marginal probability of a single item $i\in\mathcal{Y}$ from an eigendecomposition of $C$ in $O(D^{2})$ time. Similarly, given two items $i$ and $j$ we have

so we can compute arbitrary entries of $K$ in $O(D^{2})$ time. This allows us to compute, for example, pairwise marginals $\mathcal{P}(i,j\in\boldsymbol{Y})=K_{ii}K_{jj}-K_{ij}^{2}$. More generally, for a set $A\in\mathcal{Y}$, $|A|=k$, we need to compute $\frac{k(k+1)}{2}$ entries of $K$ to obtain $K_{A}$, and taking the determinant then yields $\mathcal{P}(A\subseteq\boldsymbol{Y})$. The process requires only $O(D^{2}k^{2}+k^{\omega})$ time; for small sets $|A|$ this is just quadratic in the dimension of $\phi$.

3.3 Sampling

Note that, when $\hat{V}$ contains eigenvectors of $C$, this is (up to scale) the relationship established by Proposition 3.1 between eigenvectors $\boldsymbol{\hat{v}}$ of $C$ and eigenvectors $\boldsymbol{v}$ of $L$.

The mapping in Equation (123) has several useful properties. If $\boldsymbol{v}_{1}=B^{\top}\boldsymbol{\hat{v}}_{1}$ and $\boldsymbol{v}_{2}=B^{\top}\boldsymbol{\hat{v}}_{2}$, then $\boldsymbol{v}_{1}+\boldsymbol{v}_{2}=B^{\top}(\boldsymbol{\hat{v}}_{1}+\boldsymbol{\hat{v}}_{2})$, and likewise for any arbitrary linear combination. In other words, we can perform implicit scaling and addition of the vectors in $V$ using only their preimages in $\hat{V}$. Additionally, we have

so we can compute dot products of vectors in $V$ in $O(D^{2})$ time. This means that, for instance, we can implicitly normalize $\boldsymbol{v}=B^{\top}\boldsymbol{\hat{v}}$ by updating $\boldsymbol{\hat{v}}\leftarrow\frac{\boldsymbol{\hat{v}}}{\sqrt{\boldsymbol{\hat{v}}^{\top}C\boldsymbol{\hat{v}}}}$.

We now show how these operations allow us to efficiently implement key parts of the sampling algorithm. Because the nonzero eigenvalues of $L$ and $C$ are equal, the first loop of the algorithm, where we choose in index set $J$, remains unchanged. Rather than using $J$ to construct orthonormal $V$ directly, however, we instead build the set $\hat{V}$ by adding $\frac{\boldsymbol{\hat{v}}_{n}}{\sqrt{\boldsymbol{\hat{v}}_{n}^{\top}C\boldsymbol{\hat{v}}_{n}}}$ for every $n\in J$.

In the last phase of the loop, we need to find an orthonormal basis $V_{\bot}$ for the subspace of $V$ orthogonal to a given $\boldsymbol{e}_{i}$. This requires two steps. In the first, we subtract a multiple of one of the vectors in $V$ from all of the other vectors so that they are zero in the $i$th component, leaving us with a set of vectors spanning the subspace of $V$ orthogonal to $\boldsymbol{e}_{i}$. In order to do this we must be able to compute the $i$th component of a vector $\boldsymbol{v}\in V$; since $\boldsymbol{v}=B^{\top}\boldsymbol{\hat{v}}$, this is easily done by computing the $i$th column of $B$, and then taking the dot product with $\boldsymbol{\hat{v}}$. This takes only $O(D)$ time. In the second step, we use the Gram-Schmidt process to convert the resulting vectors into an orthonormal set. This requires a series of dot products, sums, and scalings of vectors in $V$; however, as previously argued all of these operations can be performed implicitly. Therefore the mapping in Equation (123) allows us to implement the final line of the second loop using only tractable computations on vectors in $\hat{V}$.

All that remains, then, is to efficiently choose an item $i$ according to the distribution

in the first line of the while loop. Simplifying, we have

Thus the required distribution can be computed in time $O(NDk)$, where $k=|\hat{V}|$. The complete dual sampling algorithm is given in Algorithm 3; the overall runtime is $O(NDk^{2}+D^{2}k^{3})$.

4 Random projections

As we have seen, dual DPPs allow us to deal with settings where $N$ is too large to work efficiently with $L$ by shifting the computational focus to the dual kernel $C$, which is only $D\times D$. This is an effective approach when $D\ll N$. Of course, in some cases $D$ might also be unmanageably large, for instance when the diversity features are given by word counts in natural language settings, or high-resolution image features in vision.

To address this problem, we describe a method for reducing the dimension of the diversity features while maintaining a close approximation to the original DPP model. Our approach is based on applying random projections, an extremely simple technique that nonetheless provides an array of theoretical guarantees, particularly with respect to preserving distances between points (Vempala, 2004). A classic result of Johnson and Lindenstrauss (1984), for instance, shows that high-dimensional points can be randomly projected onto a logarithmic number of dimensions while approximately preserving the distances between them. More recently, Magen and Zouzias (2008) extended this idea to the preservation of volumes spanned by sets of points. Here, we apply the connection between DPPs and spanned volumes to show that random projections allow us to reduce the dimensionality of $\phi$, dramatically speeding up inference, while maintaining a provably close approximation to the original, high-dimensional model. We begin by stating a variant of Magen and Zouzias’ result.

(Adapted from Magen and Zouzias (2008)) Let $X$ be a $D\times N$ matrix. Fix $k<N$ and $0<\epsilon,\delta<1/2$, and set the projection dimension

Let $G$ be a $d\times D$ random projection matrix whose entries are independently sampled from $\mathcal{N}(0,\frac{1}{d})$, and let $X_{Y}$, where $Y\subseteq\{1,2,\dots,N\}$, denote the $D\times|Y|$ matrix formed by taking the columns of $X$ corresponding to indices in $Y$. Then with probability at least $1-\delta$ we have, for all $Y$ with cardinality at most $k$:

Lemma 3.2 says that, with high probability, randomly projecting to

dimensions is sufficient to approximately preserve all volumes spanned by $k$ columns of $X$. We can use this result to bound the effectiveness of random projections for DPPs.

In order to obtain a result that is independent of $D$, we will restrict ourselves to the portion of the distribution pertaining to subsets $Y$ with cardinality at most a constant $k$. This restriction is intuitively reasonable for any application where we use DPPs to model sets of relatively small size compared to $N$, which is common in practice. However, formally it may seem a bit strange, since it implies conditioning the DPP on cardinality. In Section 5 we will show that this kind of conditioning is actually very practical and efficient, and Theorem 3.3, which we prove shortly, will apply directly to the $k$-DPPs of Section 5 without any additional work.

For now, we will seek to approximate the distribution $\mathcal{P}^{\leq k}(Y)=\mathcal{P}(\boldsymbol{Y}=Y\ |\ |\boldsymbol{Y}|\leq k)$, which is simply the original DPP conditioned on the cardinality of the modeled subset:

where $\phi(Y)$ denotes the $D\times|Y|$ matrix formed from columns $\phi_{i}$ for $i\in Y$. Our main result follows.

Let $\mathcal{P}^{\leq k}(Y)$ be the cardinality-conditioned DPP distribution defined by quality model $q$ and $D$-dimensional diversity feature function $\phi$, and let

be the cardinality-conditioned DPP distribution obtained by projecting $\phi$ with $G$. Then for projection dimension $d$ as in Equation (129), we have

with probability at least $1-\delta$. Note that $e^{6k\epsilon}-1\approx 6k\epsilon$ when $k\epsilon$ is small.

The theorem says that for $d$ logarithmic in $N$ and linear in $k$, the $L_{1}$ variational distance between the original DPP and the randomly projected version is bounded. In order to use Lemma 3.2, which bounds volumes of parallelepipeds, to prove this bound on determinants, we will make use of the following relationship:

In order to handle the conditional DPP normalization constant

we also must adapt Lemma 3.2 to sums of determinants. Finally, for technical reasons we will change the symmetry of the upper and lower bounds from the sign of $\epsilon$ to the sign of the exponent. The following lemma gives the details.

Under the definitions and assumptions of Lemma 3.2, we have, with probability at least $1-\delta$,

where the first inequality holds with probability at least $1-\delta$ by Lemma 3.2, and the third follows from the fact that $(1-\epsilon)(1+2\epsilon)\geq 1$ (since $\epsilon<1/2$), thus $(1-\epsilon)^{2k}\geq(1+2\epsilon)^{-2k}$. The upper bound follows directly:

Recall the matrix $B$, whose columns are given by $B_{i}=q_{i}\phi_{i}$. We have

where the first inequality follows from Lemma 3.2 and Lemma 137, which hold simultaneously with probability at least $1-\delta$, and the second follows from $(1+a)^{b}\leq e^{ab}$ for $a,b\geq 0$. ∎

By combining the dual representation with random projections, we can deal simultaneously with very large $N$ and very large $D$. In fact, in Section 6 we will show that $N$ can even be exponentially large if certain structural assumptions are met. These techniques vastly expand the range of problems to which DPPs can be practically applied.

5 Alternative likelihood formulas

Recall that, in an L-ensemble DPP, the likelihood of a particular set $Y\subseteq\mathcal{Y}$ is given by

This expression has some nice intuitive properties in terms of volumes, and, ignoring the normalization in the denominator, takes a simple and concise form. However, as a ratio of determinants on matrices of differing dimension, it may not always be analytically convenient. Minors can be difficult to reason about directly, and ratios complicate calculations like derivatives. Moreover, we might want the likelihood in terms of the marginal kernel $K=L(L+I)^{-1}=I-(L+I)^{-1}$, but simply plugging in these identities yields a expression that is somewhat unwieldy.