Scalable Learning and MAP Inference for Nonsymmetric Determinantal Point Processes

1 Introduction

Determinantal point processes (DPPs) have proven useful in numerous machine learning applications. For example, recent uses include summarization Sharghi et al. (2018), recommender systems Wilhelm et al. (2018)

, neural network compression

Mariet and Sra (2016), kernel approximation Li et al. (2016), multi-modal output generation Elfeki et al. (2019), and batch selection, both for stochastic optimization Zhang et al. (2017)

and for active learning

Bıyık et al. (2019). DPPs have been applied in these cases because for each there existed a subset-selection problem, and a good solution to this problem was to select items that were high-quality but also diverse; DPPs provide a means of trading off quality with diversity in a principled way Kulesza et al. (2012).

For subset selection problems where the ground set of items to select from has cardinality $M$ , the typical DPP is parameterized by an $M \times M$ kernel matrix. Most prior work has been concerned with symmetric DPPs, where the kernel must equal its transpose. However, recent work has considered the more general class of nonsymmetric DPPs (NDPPs) and shown that these have additional useful modeling power Brunel (2018); Gartrell et al. (2019). In particular, NDPPs allow modeling of positive correlations between items, where the presence of item $i$

in the selected set increases the probability that some other item

j

will also be selected. Symmetric DPPs cannot capture such positive correlations. There are many intuitive examples of how positive correlations can be of practical importance. For example, consider a product recommendation task for a website, where a user has a camera in her shopping cart, and the goal is to display several other items that she might purchase. Relative to an empty cart, the presence of the camera probably increases the probability that the person will buy an accessory like a tripod. Although NDPPs can model such behavior, scalable approaches for NDPP learning and MAP inference have not been studied.

Here, we build on prior NDPP work and make the following contributions:

Learning: We propose a decomposition for NDPP kernels that reduces the complexity of learning from cubic Gartrell et al. (2019) to linear in $M$ , the size of the DPP ground set.

MAP inference: For the problem of finding the highest-probability subset under a DPP (the “MAP inference” problem), we analyze a standard greedy subset selection algorithm and show that, for low-rank NDPPs, it can be run in time linear in the size of the DPP ground set. We also describe how the same algorithm can be used to do conditioning for next-item prediction tasks. We then prove a lower bound on greedy MAP performance for all NDPPs.

We combine the above contributions through experiments that involve learning NDPP kernels and applying MAP inference to these kernels for several real-world datasets. These experiments demonstrate that our scalable decomposition matches or exceeds the predictive performance of prior work. We also show that, for small synthetic kernels generated from learned kernels, the proposed greedy MAP algorithm has an approximation quality that often substantially exceeds the lower bound, with near-optimal performance in many cases.

2 Background

Consider a finite set $Y = {1, 2, \dots, M}$ of cardinality $M$ , which we will also denote by $[[M]]$ . A DPP on $[[M]]$

defines a probability distribution over all its

2^{M}

subsets. It is parameterized by a matrix

L \in R^{M \times M}

, called the kernel, such that the probability of each subset

Y \subseteq [[M]]

is proportional to the determinant of its corresponding principal submatrix:

Pr (Y) \propto det (L_{Y})

. The normalization constant for this distribution can be expressed as a single

M \times M

determinant:

\sum_{Y \subseteq [[M]]} det (L_{Y}) = det (L + I)

(Kulesza et al., 2012, Theorem 2.1). Hence,

Pr (Y) = det (L_{Y}) / det (L + I)

. We will use

P_{L}

to denote this distribution.

For intuition about the kernel parameters, notice that the probabilities of singletons ${i}$ and ${j}$ are proportional to $L_{i i}$ and $L_{j j}$ , respectively. Hence, it is common to think of $L$ ’s diagonal as representing item qualities. The probability of the set ${i, j}$ is proportional to $det (L_{{i, j}}) = L_{i i} L_{j j} - L_{i j} L_{j i}$ . Thus, if $L_{i j} L_{j i} < 0$ , then $i$ and $j$ are positively correlated. Similarly, if $L_{i j} L_{j i} > 0$ , then $i$ and $j$ are negatively correlated. Therefore, off-diagonal terms determine item correlations.

In order to ensure that $P_{L}$ defines a probability distribution, all principal minors of $L$ must be non-negative: $det (L_{J}) \geq 0$ . Matrices that satisfy this property are called $P_{0}$ -matrices (Fang, 1989, Definition 1). There is no known generation method or matrix decomposition that fully covers the space of all $P_{0}$ matrices, although there are many that partially cover the space Tsatsomeros (2004).

One common partial solution is to use a decomposition that covers the space of symmetric $P_{0}$ matrices. By restricting to the space of symmetric matrices, one can exploit the fact that $L \in P_{0}$ if and only if $L$ is positive semidefinite (PSD) Prussing (1986). (Recall that a matrix $L \in R^{M \times M}$ is defined to be PSD if and only if $x^{⊤} L x \geq 0$ , for all $x \in R^{M}$

.) Any symmetric PSD matrix can be written as the Gramian matrix of some set of vectors:

L := V V^{⊤}

, where

V \in R^{M \times K}

. Hence, the

V V^{⊤}

decomposition provides an easy means of generating the entire space of symmetric

P_{0}

matrices. It also has a nice intuitive interpretation, if we view the

i

-th row of

V

as a length-

K

feature vector describing item

i

Unfortunately, the symmetry requirement severely limits the types of correlations that a DPP can capture. A symmetric model is able to capture only nonpositive correlations between items, since $L_{i j} L_{j i} = L_{i j}^{2} \geq 0$ , whereas a nonsymmetric $L$ can capture both negative and positive correlations (see (Gartrell et al., 2019, Section 2.1) for more intuition). To expand coverage to nonsymmetric matrices in $P_{0}$ , it’s natural to consider nonsymmetric PSD matrices. In what follows, we denote by $P_{0}^{+}$ the set of all nonsymmetric (and symmetric) PSD matrices. Any nonsymmetric PSD matrix is in $P_{0}$ (Gartrell et al., 2019, Lemma 1), so $P_{0}^{+} \subseteq P_{0}$ . However, unlike in the symmetric case, the set of nonsymmetric PSD matrices does not fully cover the set of nonsymmetric $P_{0}$ matrices. For example, consider

L = (\begin{matrix} 1 & 5 / 3 1 / 2 & 1 \end{matrix}) with det (L_{{1}}), det (L_{{2}}), det (L_{{1, 2}}) \geq 0, but x^{⊤} L x < 0 for x = (\begin{matrix} - 1 1 \end{matrix}) .

Still, nonsymmetric PSD matrices cover a large enough portion of the $P_{0}$ space to be useful in practice, as evidenced by the experiments of Gartrell et al. (2019). This work covered the $P_{0}^{+}$ space by using the following decomposition: $L := S + A$ , with $S := V V^{⊤}$ for $V \in R^{M \times K}$ , and $A := B C^{⊤} - C B^{⊤}$ for $B, C \in R^{M \times K}$ . This decomposition makes use of the fact that any matrix $L$ can be decomposed uniquely as the sum of a symmetric matrix $S$

and a skew-symmetric matrix

A

. (To see this, let

S = (L + L^{⊤}) / 2

and

A = (L - L^{⊤}) / 2

.) All skew-symmetric matrices

A

are trivially PSD, since

x^{⊤} A x = 0

for all

x \in R^{M}

. Hence, the

L

here is guaranteed to be PSD simply because its

S

uses the standard Gramian decomposition

V V^{⊤}

In this work we will also only consider $P_{0}^{+}$ , and leave to future work the problem of finding tractable ways to cover the rest of $P_{0}$ . We propose a new decomposition of $L$ that allows for more scalable learning. As in prior DPP work, our decomposition has inner dimension $K$ that could be as large as $M$ , but is usually much smaller in practice. Our algorithms work well for modest values of $K$ . In cases where the natural $K$

is larger (e.g., natural language processing applications where

K

might be the number of words), random projections can often be used to significantly reduce

K

Gillenwater et al. (2012a).

3 New kernel decomposition and scalable learning

Prior work on NDPPs proposed a maximum likelihood estimation (MLE) learning algorithm

Gartrell et al. (2019). Due to that work’s particular kernel decomposition, this learning approach has complexity cubic in the number of items

M

. Here, we propose an alternative decomposition that reduces the complexity to linear in

M

We begin by showing that our new decomposition covers the space of $P_{0}^{+}$ matrices. Before diving in, let us define $Σ_{i} := (\begin{matrix} 0 & λ_{i} - λ_{i} & 0 \end{matrix})$ as shorthand for a $2 \times 2$ block matrix with zeros on-diagonal and opposite values off-diagonal. Then, our proposed decomposition is as follows:

L := S + A, with S := V V^{⊤} and A := B C B^{⊤}

(1)

where $V, B \in R^{M \times K}$ , and $C \in R^{K \times K}$ is a block-diagonal matrix with some diagonal blocks of the form $Σ_{i}$ , with $λ_{i} > 0$ , and zeros elsewhere. The following lemma shows that this decomposition covers the space of $P_{0}^{+}$ matrices.

Lemma 1.

For $k \in [[M]]$ , let $ℓ \leq M$ be an even integer and let $A \in R^{M \times M}$ be a skew-symmetric matrix with rank $ℓ$ . Then, there exist $B \in R^{M \times ℓ}$ and positive numbers $λ_{1}, \dots, λ_{ℓ / 2}$ , such that $A = B C B^{⊤}$ , where $C \in R^{ℓ \times ℓ}$ is the block-diagonal matrix with $(ℓ / 2)$ diagonal blocks of size 2 given by $Σ_{i}$ , $i = 1, \dots, ℓ / 2$ .

The proof of lemma:skew-sym-decomp and all subsequent results can be found in sec:proof. With this decomposition in hand, we now proceed to show that it can be used for linear-time MLE learning. To do so, we must show that corresponding NDPP log-likelihood objective and gradient can be computed in time linear in $M$ . Given a collection of $n$ observed subsets ${Y_{1}, . . ., Y_{n}}$ composed of items from $Y = [[M]]$ , the full formulation of the regularized log-likelihood is:

ϕ (V, B, C)

= n \sum i = 1 log det (V_{Y_{i}} V_{Y_{i}}^{⊤} + B_{Y_{i}} C B_{Y_{i}}^{⊤}) - log det (V V^{⊤} + B C B^{⊤} + I) - R (V, B),

(2)

where $V_{Y_{i}} \in R^{| Y_{i} | \times K}$ denotes a matrix composed of the rows of $V$ that correspond to the items in $Y_{i}$ . The regularization term, $R (V, B)$ , is defined as follows:

R (V, B)

= α M \sum i = 1 \frac{1}{λ_{i}} ∥ v_{i} ∥_{2}^{2} + β M \sum i = 1 \frac{1}{λ_{i}} ∥ b_{i} ∥_{2}^{2},

(3)

where $λ_{i}$ counts the number of occurrences of item $i$ in the training set, $v_{i}$ and $b_{i}$ are rows of $V$ and $B$ , respectively, and $α, β, > 0$

are tunable hyperparameters. This regularization is similar to that of prior work

Gartrell et al. (2017, 2019). We omit regularization for

C

since regularization based on item counts (popularity) cannot be directly applied to

C

. Furthermore, we observe from our experiments that regularization on

C

does not seem to be needed in practice.

theorem:loglik-complexity,theorem:grad-complexity show that computing the regularized log-likelihood and its gradient both have time complexity linear in $M$ . The complexities also depend on $K$ , the rank of the NDPP, and $K^{'}$ , the size of the largest observed subset in the data. In practice though, usually $K ≪ M$ and $K = K^{'}$ . Hence, linearity in $M$ means that we can efficiently perform learning for datasets with very large ground sets, which is impossible with the cubic-complexity $L$ decomposition in prior work Gartrell et al. (2019). We briefly note here that a key component in the proof of theorem:loglik-complexity is use of the Woodbury matrix identity. This allows us to transform an $O (M^{3})$ operation into an $O (M K^{2})$ one, by replacing the computation of a size- $M$ $det$ with $O (M K^{2})$ matrix multiplications and a size- $K$ $det$ . A similar approach was also recently adopted in work on normalizing flows for deep generative models Lu and Huang (2020).

Theorem 1.

Given a nonsymmetric low-rank DPP parameterized by $V$ of rank $K$ , $B$ of rank $K$ , and a $K \times K$ matrix $C$ , where $L = V V^{⊤} + B C B^{⊤}$ , we can compute the regularized log-likelihood, eq:nonsymm-full-log-likelihood, in $O (M K^{2} + K^{3} + n K^{' 3})$ time, where $K^{'}$ is the size of the largest of the $n$ training subsets.

Theorem 2.

Given a nonsymmetric low-rank DPP parameterized by $V$ of rank $K$ , $B$ of rank $K$ , and a $K \times K$ matrix $C$ , where $L = V V^{⊤} + B C B^{⊤}$ , we can compute the gradient of the regularized log-likelihood in eq:nonsymm-full-log-likelihood in $O (M K^{2} + K^{3} + n K^{' 3})$ time, where $K^{'}$ is the size of the largest of the $n$ training subsets.

To further simplify learning and MAP inference, we set $B = V$ , which results in $L = V V^{⊤} + V C V^{⊤} = V (I + C) V^{⊤}$ . This change also simplifies regularization, so that we only perform regularization on $V$ , as indicated in the first term of eq:regularization, leaving us with the single regularization hyperparameter of $α$ . While setting $B = V$ restricts the class of nonsymmetric $L$ kernels that can be represented, we find in practice that this does not adversely impact prediction quality.

To compensate for the restriction imposed by setting $B = V$ , we also relax the block-diagonal structure imposed on $C$ , so that we learn a full skew-symmetric $K \times K$ matrix $C$ . We empirically observe that this relaxation of the block-diagonal structure is needed to ensure learning of nonsymmetric $L$ kernels with better predictive performance. To ensure that $C$ and thus $A$ is skew-symmetric, we set $C = D - D^{T}$ , where $D \in R^{K \times K}$ .

4 MAP inference

After learning an NDPP, one can then use it to infer the most probable item subsets in various situations. Several inference algorithms have been well-studied for symmetric DPPs, including sampling Kulesza and Taskar (2011); Anari et al. (2016); Li et al. (2016); Launay et al. (2018); Gillenwater et al. (2019); Poulson (2019); Dereziński (2019) and MAP inference Gillenwater et al. (2012b); Han et al. (2017); Chen et al. (2018); Han and Gillenwater (2020). We focus on MAP inference:

argmax Y \subseteq Y det (L_{Y}) such that | Y | = k,

(4)

for cardinality budget $k \leq M$ . MAP inference for DPPs is known to be NP-hard even in the symmetric case Ko et al. (1995); Kulesza et al. (2012). For symmetric DPPs, one usually approximates the MAP via the standard greedy algorithm for submodular maximization Nemhauser et al. (1978). First, we describe how to efficiently implement this for NDPPs. Then, in sec:map-approx-guarantee we prove a lower bound on its approximation quality. To the best of our knowledge, this is the first investigation of how to apply the greedy algorithm to NDPPs.

Greedy begins with an empty set and repeatedly adds the item that maximizes the marginal gain until the chosen set is size $k$ . Here, we aim to design an efficient greedy algorithm when the NDPP kernel is given by a low-rank matrix. For generality, in what follows we write the kernel as $L = B C B^{⊤}$ , since one can easily rewrite our matrix decomposition, as well as that of Gartrell et al. (2019), to take this form. For example, for our decomposition: . Using Schur’s determinant identity, we first observe that, for $Y \subseteq [[M]]$ and $i \in [[M]]$ , the marginal gain of a NDPP can be written as

\frac{det (L_{Y \cup {i}})}{det (L_{Y})}

= L_{i i} - L_{i Y} (L_{Y})^{- 1} L_{Y i} = b_{i} C b_{i}^{⊤} - b_{i} C (B_{Y}^{⊤} (B_{Y} C B_{Y}^{⊤})^{- 1} B_{Y}) C b_{i}^{⊤},

(5)

where $b_{i} \in R^{1 \times K}$ and $B_{Y} \in R^{| Y | \times K}$ . A naïve computation of eq:bilinear-marginalgain is $O (K^{2} + k^{3})$ , since we must invert a $| Y | \times | Y |$ matrix, where $| Y | \leq k$ . However, one can compute eq:bilinear-marginalgain more efficiently by observing that it actually can be expressed as a rank- $| Y |$ matrix and hence computed in $O (K^{2})$ time.

Lemma 2.

Given $B \in R^{M \times K}$ , $C \in R^{K \times K}$ , and $Y = {a_{1}, \dots, a_{k}} \subseteq [[M]]$ , let $b_{i} \in R^{1 \times K}$ be the $i$ -th row in $B$ and $B_{Y} \in R^{| Y | \times K}$ be a matrix containing rows in $B$ indexed by $Y$ . Then, it holds that

B_{Y}^{⊤} (B_{Y} C B_{Y}^{⊤})^{- 1} B_{Y} = k \sum j = 1 p_{j}^{⊤} q_{j},

(6)

where row vectors $p_{j}, q_{j} \in R^{1 \times K}$ for $j = 1, \dots, k$ satisfy $p_{1} = b_{a_{1}} / (b_{a_{1}} C b_{a_{1}}^{⊤})$ , $q_{1} = b_{a_{1}}$ , and

p_{j + 1} = \frac{b_{a_{j}} - b_{a_{j}} C^{⊤} \sum_{i = 1}^{j} q_{i}^{⊤} p_{i}}{b_{a_{j}} C (b_{a_{j}} - b_{a_{j}} C^{⊤} \sum_{i = 1}^{j} q_{i}^{⊤} p_{i})^{⊤}}, q_{j + 1} = b_{a_{j}} - b_{a_{j}} C j \sum i = 1 p_{i}^{⊤} q_{i} .

(7)

Plugging eq:bilinear-inverse into eq:bilinear-marginalgain, the marginal gain with respect to $Y \cup {a}$ can be computed by simply updating from the previous gain with respect to $Y$ . That is,

	$\frac{det (L_{Y \cup {a, i}})}{det (L_{Y \cup {a}})}$	$= b_{i} C b_{i}^{⊤} - \| Y \| + 1 \sum j = 1 (b_{i} C p_{j}^{⊤}) (b_{i} C^{⊤} q_{j}^{⊤})$		(8)
		$= \frac{det (L_{Y \cup {i}})}{det (L_{Y})} - (b_{i} C p_{\| Y \| + 1}^{⊤}) (b_{i} C^{⊤} q_{\| Y \| + 1}^{⊤}) .$		(9)

The marginal gains when $Y = \emptyset$ are equal to diagonals of $L$ and require $O (M K^{2})$ operations. Then, computing the update terms in eq:marginalgainupdate for all $i \in [[M]]$ needs $O (M K)$ operations. Since the total number of updates is $k$ , the overall complexity becomes $O (M K^{2} + M K k)$ . We provide a full description of the implied greedy algorithm for low-rank NDPPs in alg:greedylowrank.

1:Input:

B \in R^{M \times K}

C \in R^{K \times K}

, the cardinality

k

▹

And

{a_{1}, \dots, a_{k}}

for conditioning

2:initialize

P \leftarrow []

Q \leftarrow []

and

Y \leftarrow \emptyset

Δ_{i} \leftarrow b_{i} C b_{i}^{⊤}

for

i \in [[M]]

where

b_{i} \in R^{1 \times K}

is the

i

-th row in

B

a \leftarrow {argmax}_{i} Δ_{i}

and

Y \leftarrow Y \cup {a}

▹

a \leftarrow a_{1}

for conditioning

5:while

| Y | \leq k

p \leftarrow (b_{a} - b_{a} C^{⊤} Q^{⊤} P) / Δ_{a}

q \leftarrow b_{a} - b_{a} C P^{⊤} Q

P \leftarrow [P; p]

and

Q \leftarrow [Q; q]

9: for

i \in [[M]]

10:

a \leftarrow {argmax}_{i} Δ_{i}

and

Y \leftarrow Y \cup {a}

▹

a \leftarrow a_{| Y | + 1}

for conditioning

11:end while

12:return

Y

▹

return

{Δ_{i}}_{i = 1}^{M}

for conditioning

Algorithm 1 Greedy MAP inference for low-rank NDPPs

Low-rank DPP Models	MLE Learning	MAP Inference
Symmetric DPP Gartrell et al. (2017)	$O (M K^{2} + n K^{3})$	$O (M K k + M K^{2})$
Nonsymmetric DPP Gartrell et al. (2019)	$O (M^{3} + M K^{2} + n K^{3})$	$O (M K k + M K^{2})$
Scalable nonsymmetric DPP (this work)	$O (M K^{2} + n K^{3})$	$O (M K k + M K^{2})$

Table 1: Runtime complexities for several DPP models. Our model and the symmetric DPP model Gartrell et al. (2017) can perform both tasks in time linear in the size of ground set

M

, but ours is a more general model that can capture positive as well as negative item correlations.

tab:runtimes summarizes the runtime complexitiy of our methods and those of previous work. We also note that memory required for MAP inference is $O (M K)$ in all cases. For learning, the memory required is $O (M^{2})$ for the NDPPs from prior work Gartrell et al. (2019), but only $O (M K)$ for low-rank symmetric DPPs Gartrell et al. (2017) and our proposed scalable NDPPs.

4.1 Approximation guarantee for greedy NDPP MAP inference

As mentioned above, alg:greedylowrank is an instantiation of the standard greedy algorithm used for submodular maximization Nemhauser et al. (1978). This algorithm has a $(1 - 1 / e)$ -approximation guarantee for the problem of maximizing nonnegative, monotone submodular functions. While the function $f (Y) = log det (L_{Y})$ is submodular for a symmetric PSD $L$ Kelmans and Kimelfeld (1983), it is not monotone. Often, as in Han and Gillenwater (2020)

, it is assumed that the smallest eigenvalue of

L

is greater than

1

, which guarantees montonicity. There is no particular evidence that this assumption is true for typical practical models, but nevertheless the greedy algorithm tends to perform well in practice for symmetric DPPs. Here, we prove a similar approximation guarantee that covers NDPPs as well, even though the function

f (Y) = log det (L_{Y})

is non-submodular when

L

is nonsymmetric. In sec:mapinference, we further observe that, as for symmetric DPPs, the greedy algorithm seems to work well in practice for NDPPs.

Recently, Bian et al. (2017) proposed an extension of greedy algorithm guarantees to non-submodular functions. Their result is based on the submodularity ratio and curvature of the objective function, which measure to what extent the objective has submodular properties. Leveraging this result, in thm:greedy we provide an approximation ratio for greedy MAP inference of NDPPs.

Theorem 3.

Consider a nonsymmetric low-rank DPP $L = V V^{⊤} + B C B^{⊤}$ , where $V, B$ are of rank $K$ , and $C \in R^{K \times K}$ . Given a cardinality budget $k$ , let $σ_{min}$ and $σ_{max}$

denote the smallest and largest singular values of

L_{Y}

for all

Y \subseteq [[M]]

and

| Y | \leq 2 k

. Assume that

σ_{min} > 1

. Then,

log det (L_{Y^{G}})

\geq \frac{4 (1 - e^{- 1 / 4})}{2 (log σ_{max} / log σ_{min}) - 1} log det (L_{Y^{*}})

(10)

where $Y^{G}$ is the output of alg:greedylowrank and $Y^{*}$ is the optimal solution of MAP inference in eq:mapinference.

Thus, when the kernel has a small value of $log σ_{max} / log σ_{min}$ , the greedy algorithm finds a near-optimal solution. In practice, we observe that the greedy algorithm finds a near-optimal solution even for large values of this ratio (see sec:mapinference). As remarked above, there is no evidence that the condition $σ_{min} > 1$ is usually true in practice. While this condition can be achieved by multiplying $L$ by a constant, this leads to a (potentially large) additive term in eq:greedyratio. We provide cor:greedy in sec:corr-greedy, which excludes the $σ_{min} > 1$ assumption, and quantifies this additive term.

4.2 Greedy conditioning for next-item prediction

We briefly describe here a small modification to the greedy algorithm that is necessary if one wants to use it as a tool for next-item prediction. Given a set $Y \subseteq [[M]]$ , Kulesza et al. (2012) showed that a DPP with $L$ conditioned on the inclusion of the items in the set $Y$ forms another DPP with kernel $L^{Y} := L_{¯ Y} - L_{¯ Y, Y} L_{Y}^{- 1} L_{¯ Y, Y}$ where $¯ Y = [[M]] ∖ Y$ . The singleton probability $Pr (Y \cup {i} ∣ Y) \propto L_{i i}^{Y}$ can be useful for doing next-item prediction. We can use the same machinery from the greedy algorithm’s marginal gain computations to effectively compute these singletons. More concretely, suppose that we are doing next-item prediction as a shopper adds items to a digital cart. We predict the item that maximizes the marginal gain, conditioned on the current cart contents (the set $Y$ in the greedy algorithm). When the shopper adds the next item to their cart, we update $Y$ to include this item, rather than our predicted item (line 10 in alg:greedylowrank). We then iterate until the shopper checks out. The comments on the righthand side of alg:greedylowrank summarize this procedure. The runtime of this prediction is the same that of the greedy algorithm, $O (M K^{2} + M K | Y |)$ . We note that this cost is comparable to that of an approach based on the DPP dual kernel from prior work Mariet et al. (2019), which has $O (M K^{2} + K^{3} + | Y |^{3})$ complexity. However, since it is non-trivial to define the dual kernel for NDPPs, the greedy algorithm may be the simpler choice for next-item prediction for NDPPs.

5 Experiments

Code for all experiments is available at
https://github.com/cgartrel/nonsymmetric-DPP-learning/tree/scalable.

5.1 Datasets

We perform experiments on several real-world public datasets composed of online shopping baskets:

Amazon Baby Registries: This dataset consists of registries or "baskets" of baby products, and has been used in prior work on DPP learning Gartrell et al. (2016, 2019); Gillenwater et al. (2014); Mariet and Sra (2015). The registries contain items from 15 different categories, such as “apparel”, with a catalog of up to 100 items per category. Our evaluation mirrors that of Gartrell et al. (2019); we evaluate on the popular apparel category, which contains 14,970 registries, as well as on a dataset composed of the three most popular categories: apparel, diaper, and feeding, which contains a total of 31,218 registries.
UK Retail: This dataset Chen et al. (2012) contains baskets representing transactions from an online retail company that sells unique all-occasion gifts. We omit baskets with more than 100 items, leaving us with a dataset containing 19,762 baskets drawn from a catalog of 3,941 products. Baskets containing more than 100 items are in the long tail of the basket-size distribution of the data, so omitting larger baskets is reasonable, and allows us to use a low-rank factorization of the DPP with $K = 100$ .
Instacart: This dataset Instacart (2017) contains baskets purchased by Instacart users. Again, we omit baskets with more than 100 items, leaving us with 3.2 million baskets and a catalog of 49,677 products.

5.2 Experimental setup and metrics

We use a small held-out validation set, consisting of 100 randomly-selected baskets, for tracking convergence during training and for tuning hyperparameters. A random selection of 2000 of the remaining baskets are used for testing, and the rest are used for training. Convergence is reached during training when the relative change in validation log-likelihood is below a predetermined threshold. We use PyTorch with Adam

Kingma and Ba (2015) for optimization.

Subset expansion task. We use greedy conditioning to do next-item prediction, as described in subsec:greedy-conditioning. We compare methods using a standard recommender system metric: mean percentile rank (MPR) Hu et al. (2008); Li et al. (2010). MPR of 50 is equivalent to random selection; MPR of 100 means that the model perfectly predicts the next item. See sec:MPR for a complete description of the MPR metric.

Subset discrimination task. We also test the ability of a model to discriminate observed subsets from randomly generated ones. For each subset in the test set, we generate a subset of the same length by drawing items uniformly at random (and we ensure that the same item is not drawn more than once for a subset). We compute the AUC for the model on these observed and random subsets, where the score for each subset is the log-likelihood that the model assigns the subset.

Amazon: Apparel (

M = 100

)

Amazon: 3-category (

M = 300

)

Metric

Sym

Nonsym

Scalable nonsym

Sym

Nonsym

Scalable nonsym

MPR

77.42

\pm

1.12

80.32

\pm

0.75

84.86

\pm

1.51

60.61

\pm

0.94

75.09

\pm

0.85

75.08

\pm

1.64

AUC

0.66

\pm

0.01

0.73

\pm

0.01

0.71

\pm

0.03

0.70

\pm

0.01

0.79

\pm

0.01

0.80

\pm

0.04

test log-likelihood

-10.31

-9.66

-9.64

-18.11

-17.02

-17.22

UK Retail (

M = 3,941

)

Instacart (

M = 49,677

)

Metric

Sym

Nonsym

Scalable nonsym

Sym

Nonsym

Scalable nonsym

MPR

76.79

\pm

0.60

79.45

\pm

0.57

81.64

\pm

1.22

91.01

\pm

1.22

92.90

\pm

0.49

AUC

0.57

\pm

0.001

0.65

\pm

0.01

0.76

\pm

0.02

0.83

\pm

0.01

0.84

\pm

0.01

test log-likelihood

-120.47

-108.67

-106

-70.81

-69.12

Table 2: Average MPR, AUC, and test log-likelihood for all datasets, for the low-rank symmetric DPP Gartrell et al. (2017), low-rank NDPP Gartrell et al. (2019)

, and our scalable NDPP models. MPR and AUC results show 95% confidence estimates obtained via bootstrapping. Bold values indicate improvement over the symmetric low-rank DPP outside of the confidence interval. See sec:hyperparams for the hyperparameter settings used in these experiments. The baseline NDPP model cannot be feasibly trained on Instacart, as memory and computational costs are prohibitive due to its large ground set size.

5.3 Predictive performance results for learning

Since the focus of our work is on improving NDPP scalability, we use the low-rank symmetric DPP Gartrell et al. (2017) and the low-rank NDPP of prior work Gartrell et al. (2019) as baselines for our experiments. tab:predictive-qual compares these DPP approaches and our scalable low-rank NDPP. We see that our scalable NDPP matches or exceeds the predictive quality of the baseline NDPP. Notice that our scalable NDPP opens to the door to training on datasets with large $M$ , such as the Instacart dataset, which is infeasible for the baseline NDPP due to very high memory and compute costs. As expected, we also empirically observe that the scalable NDPP trains far faster than the NDPP for datasets with large ground sets. For example, the per-iteration gradient update of scalable NDPP is $8 \times$ faster than that of the decomposition of Gartrell et al. (2019) on the UK dataset. See sec:training-time for a comparison of overall training times.

5.4 Performance results for MAP inference

We run various approximatation algorithms for MAP inference, including the greedy algorithm (alg:greedylowrank), stochastic greedy algorithm Mirzasoleiman et al. (2015), MCMC-based DPP sampling Li et al. (2016), and greedy local search Kathuria and Deshpande (2016). The stochastic greedy algorithm computes marginal gains of a few items chosen uniformly at random and selects the best among them. The MCMC sampling begins with a random subset $Y$ of size $k$ and picks $i \in Y$ and $j \notin Y$ uniformly at random. Then, it swaps them with probability $min (1, det (L_{Y \cup {i} ∖ {j}}) / det (L_{Y}))$ and iterates this process. The greedy local search algorithm Kathuria and Deshpande (2016) starts from the output from the greedy algorithm, $Y^{G}$ , and replaces $i \in Y^{G}$ with $j \notin Y^{G}$ that gives the maximum improvement, if such $i, j$ exist. This replacement process iterates until no improvement exists, or at most $k^{2} log (10 k)$ steps have been completed, to guarantee a tight approximation Kathuria and Deshpande (2016). We use greedy local search as a baseline since it always returns a better solution than greedy. However, it is the slowest among all algorithms, as its time complexity is $O (M K k^{5} log k)$ . We choose $k = 10$ , and provide more details of all algorithms in sec:benchmark.

To evaluate the performance of MAP inference, we report the relative log-determinant ratio, defined as $∣ ∣ \frac{log det (L_{Y^{*}}) - log det (L_{Y})}{log det (L_{Y^{*}})} ∣ ∣$ where $Y$ is the output of benchmark algorithms and $Y^{*}$ is the greedy local search result. Results are reported in tab:mapinference. We observe that the greedy algorithm achieves performance close to that of the significantly more expensive greedy local search algorithm, with relative errors of up to $0.061$ . Stochastic greedy and MCMC sampling have significantly larger errors.

Algorithms

Amazon: Apparel

Amazon: 3-category

UK Retail

Instacart

Greedy (alg:greedylowrank)

0.0196

\pm

0.0020

0.0613

\pm

0.0027

0.0498

\pm

0.0017

0.0085

\pm

0.0005

Stochastic greedy Mirzasoleiman et al. (2015)

0.1296

\pm

0.0041

0.1716

\pm

0.0041

0.1526

\pm

0.0028

0.1408

\pm

0.0044

MCMC sampling Li et al. (2016)

0.5437

\pm

0.0083

0.7640

\pm

0.0092

0.8930

\pm

0.0085

3.2646

\pm

0.0401

Table 3: Average relative error and 95% confidence intervals of MAP inference algorithms on NDPPs learned from real-world datasets. For all datasets, we evaluate

10

kernels learned with different initializations, and run

100

random trials for stochastic greedy and MCMC sampling. All errors are relative to greedy local search.

5.5 Performance guarantee for greedy MAP inference

The matrices learned on real datasets are too large to compute the exact MAP solution, but we can compute exact MAP for small matrices. In this section, we explore the performance of the greedy algorithm studied in thm:greedy for $5 \times 5$ synthetic kernel matrices. More formally, we first pick $K = 3$ singular values $s_{1}, s_{2}, s_{3}$ from a kernel learned for the “Amazon: 3-category” dataset (a plot of these singular values can be seen in fig:approxratio:singular) and generate $L = V_{1} d i a g ([s_{1}, s_{2}, s_{3}]) V_{2}^{⊤}$ , where $V_{1}, V_{2} \in R^{5 \times 3}$ are random orthonormal matrices. To ensure that $L$ is a $P_{0}$ matrix, we repeatedly sample $V_{1}, V_{2}$ until all principal minors of $L$ are nonnegative. We also evaluate the performance of the symmetric DPP, where the kernel matrices are generated similarly to the NDPP, except we set $V_{1} = V_{2}$ . We set $k = 3$ and generate $10,000$ random kernels for both symmetric DPPs and NDPPs.

The results for symmetric and nonsymmetric DPPs are shown in fig:approxratio:sym and fig:approxratio:nonsym, respectively. We plot the approximation ratio of alg:greedylowrank, i.e., $log det (L_{Y^{G}}) / log det (L_{Y^{*}})$ , with respect to $log (σ_{max} / σ_{min})$ , from cor:greedy. We observe that the greedy algorithm for both often shows approximation ratios close to $1$ . However, the worst-case ratio for NDPPs is worse than that of symmetric DPPs; $log det (L_{Y})$ for $L \in P_{0}^{+}$ is non-submodular, and the greedy algorithm with a nonsubmodular function does not have as tight of a worst-case bound as in the symmetric case.

6 Conclusion

We have presented a new decomposition for nonsymmetric DPP kernels that can be learned in time linear in the size of the ground set, which is a significant improvement over the cubic complexity of prior work. Empirical results indicate that this decomposition matches or exceeds the predictive performance of the prior decomposition. We have also proved a lower bound on the quality of the greedy MAP approximation algorithm for nonsymmetric DPPs, and showed how to implement it efficiently. For future work, we will investigate tightening the performance bound for the greedy MAP algorithm, work on developing intuition about the meaning of the parameters in the $C$ matrix, and consider kernel decompositions that cover other parts of the nonsymmetric $P_{0}$ space.

Broader Impact

In general, we feel that our work moves in a positive direction by decreasing the storage and computation costs of learning NDPPs. However, in terms of broader impact, concerns related to this work are similar to those of other recommender system work (Milano et al., 2020). When applying our methods to learn kernels from user data, we recommend employing a technique such as differentially-private SGD (Abadi et al., 2016) to help prevent user data leaks, and adjusting the weights on training examples to balance the impact of sub-groups of users so as to make the final kernel as fair as possible.

References

[1] M. Abadi, A. Chu, I. Goodfellow, H. B. McMahan, I. Mironov, K. Talwar, and L. Zhang (2016) Deep Learning with Differential Privacy. In Conference on Computer and Communications Security, Cited by: Broader Impact.
[2] N. Anari, S. O. Gharan, and A. Rezaei (2016) Monte Carlo Markov Chain Algorithms for Sampling Strongly Rayleigh Distributions and Determinantal Point Processes. In Conference on Learning Theory (COLT), Cited by: §4.
[3] A. A. Bian, J. M. Buhmann, A. Krause, and S. Tschiatschek (2017) Guarantees for Greedy Maximization of Non-submodular Functions with Applications. In International Conference on Machine Learning (ICML), Cited by: §F.5, §4.1, Theorem 4.
[4] E. Bıyık, K. Wang, N. Anari, and D. Sadigh (2019) Batch Active Learning Using Determinantal Point Processes. arXiv:1906.07975. Cited by: §1.
[5] V. Brunel (2018) Learning Signed Determinantal Point Processes through the Principal Minor Assignment Problem. In Neural Information Processing Systems (NeurIPS), Cited by: §1.
[6] D. Chen, S. L. Sain, and K. Guo (2012) Data mining for the online retail industry: A case study of RFM model-based customer segmentation using data mining. Journal of Database Marketing & Customer Strategy Management. Cited by: item 2.
[7] L. Chen, G. Zhang, and E. Zhou (2018) Fast greedy MAP inference for Determinantal Point Process to improve recommendation diversity. In Neural Information Processing Systems (NeurIPS), Cited by: §4.
[8] M. Dereziński (2019) Fast determinantal point processes via distortion-free intermediate sampling. In Conference on Learning Theory (COLT), Cited by: §4.
[9] M. Elfeki, C. Couprie, M. Riviere, and M. Elhoseiny (2019) GDPP: Learning Diverse Generations using Determinantal Point Processes. In International Conference on Machine Learning (ICML), Cited by: §1.
[10] L. Fang (1989) On the Spectra of $P$ - and $P_{0}$ -Matrices. In Linear Algebra and its Applications, Cited by: §2.
[11] M. Gartrell, V. Brunel, E. Dohmatob, and S. Krichene (2019) Learning Nonsymmetric Determinantal Point Processes. In Neural Information Processing Systems (NeurIPS), Cited by: Appendix B, Figure 2, Appendix C, §1, §1, §2, §3, §3, §3, Table 1, §4, §4, item 1, §5.3, Table 2.
[12] M. Gartrell, U. Paquet, and N. Koenigstein (2016) Bayesian low-rank determinantal point processes. In Conference on Recommender Systems (RecSys), Cited by: Appendix B, item 1.
[13] M. Gartrell, U. Paquet, and N. Koenigstein (2017) Low-Rank Factorization of Determinantal Point Processes. In
Conference on Artificial Intelligence (AAAI)
, Cited by: Appendix B, §3, Table 1, §4, §5.3, Table 2.
[14] J. Gillenwater, A. Kulesza, E. Fox, and B. Taskar (2014) Expectation-Maximization for learning Determinantal Point Processes. In Neural Information Processing Systems (NIPS), Cited by: item 1.
[15] J. Gillenwater, A. Kulesza, Z. Mariet, and S. Vassilvtiskii (2019) A Tree-Based Method for Fast Repeated Sampling of Determinantal Point Processes. In International Conference on Machine Learning (ICML), Cited by: §4.
[16] J. Gillenwater, A. Kulesza, and B. Taskar (2012) Discovering Diverse and Salient Threads in Document Collections. In Empirical Methods in Natural Language Processing (EMNLP), Cited by: §2.
[17] J. Gillenwater, A. Kulesza, and B. Taskar (2012) Near-Optimal MAP Inference for Determinantal Point Processes. In Neural Information Processing Systems (NIPS), Cited by: §4.
[18] I. Han and J. Gillenwater (2020) MAP Inference for Customized Determinantal Point Processes via Maximum Inner Product Search. In Conference on Artificial Intelligence and Statistics (AISTATS), Cited by: §4.1, §4.
[19] I. Han, P. Kambadur, K. Park, and J. Shin (2017) Faster greedy MAP inference for determinantal point processes. In International Conference on Machine Learning (ICML), Cited by: §4.
[20] Y. Hu, Y. Koren, and C. Volinsky (2008) Collaborative Filtering for Implicit Feedback Datasets. In International Conference on Data Mining (ICDM), Cited by: §5.2.
[21] Instacart (2017) The Instacart Online Grocery Shopping Dataset. Note: Accessed May 2020 External Links: Link Cited by: item 3.
[22] T. Kathuria and A. Deshpande (2016) On sampling and greedy map inference of constrained determinantal point processes. arXiv preprint arXiv:1607.01551. Cited by: Appendix D, Table 4, §5.4.
[23] A.K. Kelmans and B.N. Kimelfeld (1983) Multiplicative submodularity of a matrix’s principal minor as a function of the set of its rows and some combinatorial applications. Discrete Mathematics. Cited by: §4.1.
[24] D. P. Kingma and J. Ba (2015) Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), Cited by: §5.2.
[25] C. Ko, J. Lee, and M. Queyranne (1995) An Exact Algorithm for Maximum Entropy Sampling. Operations Research. Cited by: §4.
[26] A. Kulesza, B. Taskar, et al. (2012) Determinantal Point Processes for Machine Learning. Foundations and Trends® in Machine Learning. Cited by: §1, §2, §4.2, §4.
[27] A. Kulesza and B. Taskar (2011) Learning determinantal point processes. In Conference on Uncertainty in Artificial Intelligence (UAI), Cited by: Appendix D, §4.
[28] C. Launay, B. Galerne, and A. Desolneux (2018) Exact Sampling of Determinantal Point Processes without Eigendecomposition. arXiv preprint arXiv:1802.08429. Cited by: §4.
[29] C. Li, S. Jegelka, and S. Sra (2016) Fast DPP Sampling for Nystrom with Application to Kernel Methods. In International Conference on Machine Learning (ICML), Cited by: Appendix D, Table 4, §1, §4, §5.4, Table 3.
[30] Y. Li, J. Hu, C. Zhai, and Y. Chen (2010) Improving One-class Collaborative Filtering by Incorporating Rich User Information. In Conference on Information and Knowledge Management (CIKM), Cited by: §5.2.
[31] Y. Lu and B. Huang (2020) Woodbury Transformations for Deep Generative Flows. arXiv preprint arXiv:2002.12229. Cited by: §F.2, §3.
[32] Z. Mariet, M. Gartrell, and S. Sra (2019) Learning Determinantal Point Processes by Sampling Inferred Negatives. In Conference on Artificial Intelligence and Statistics (AISTATS), Cited by: §4.2.
[33] Z. Mariet and S. Sra (2015) Fixed-point algorithms for learning Determinantal Point Processes. In International Conference on Machine Learning (ICML), Cited by: item 1.
[34] Z. Mariet and S. Sra (2016) Diversity Networks: Neural Network Compression Using Determinantal Point Processes. In International Conference on Learning Representations (ICLR), Cited by: §1.
[35] S. Milano, M. Taddeo, and L. Floridi (2020) Recommender Systems and Their Ethical Challenges. AI & SOCIETY. Cited by: Broader Impact.
[36] B. Mirzasoleiman, A. Badanidiyuru, A. Karbasi, J. Vondrák, and A. Krause (2015) Lazier Than Lazy Greedy. In Conference on Artificial Intelligence (AAAI), Cited by: Appendix D, Table 4, §5.4, Table 3.
[37] G. Nemhauser, L. Wolsey, and M. Fisher (1978) An Analysis of Approximations for Maximizing Submodular Set Functions I. Mathematical Programming 14 (1). Cited by: §4.1, §4.
[38] J. Poulson (2019) High-performance sampling of generic Determinantal Point Processes. arXiv preprint arXiv:1905.00165. Cited by: §4.
[39] J. E. Prussing (1986) The Principal Minor Test for Semidefinite Matrices. Journal of Guidance, Control, and Dynamics. Cited by: §2.
[40] A. Sharghi, A. Borji, C. Li, T. Yang, and B. Gong (2018) Improving Sequential Determinantal Point Processes for Supervised Video Summarization. In
Proceedings of the European Conference on Computer Vision (ECCV)
, Cited by: §1.
[41] G. Thompson (1988) Normal forms for skew-symmetric matrices and Hamiltonian systems with first integrals linear in momenta. In Proceedings of the American Mathematical Society, Cited by: §F.1.
[42] R. C. Thompson (1972) Principal submatrices IX: Interlacing inequalities for singular values of submatrices. Linear Algebra and its Applications. Cited by: Theorem 5.
[43] M. J. Tsatsomeros (2004) Generating and Detecting Matrices with Positive Principal Minors. In Focus on Computational Neurobiology, Cited by: §2.
[44] M. Wilhelm, A. Ramanathan, A. Bonomo, S. Jain, E. H. Chi, and J. Gillenwater (2018) Practical Diversified Recommendations on YouTube with Determinantal Point Processes. In Conference on Information and Knowledge Management (CIKM), Cited by: §1.
[45] C. Zhang, H. Kjellström, and S. Mandt (2017) Determinantal Point Processes for Mini-Batch Diversification. In Conference on Uncertainty in Artificial Intelligence (UAI), Cited by: §1.

Appendix A Mean Percentile Rank

We begin our definition of MPR by defining percentile rank (PR). First, given a set $J$ , let $p_{i, J} = Pr (J \cup {i} ∣ J)$ . The percentile rank of an item $i$ given a set $J$ is defined as

{PR}_{i, J} = \frac{\sum_{i^{'} \notin J} 1 (p_{i, J} \geq p_{i^{'}, J})}{| Y ∖ J |} \times 100 %

where $Y ∖ J$ indicates those elements in the ground set $Y$ that are not found in $J$ .

For our evaluation, given a test set $Y$ , we select a random element $i \in Y$ and compute ${PR}_{i, Y ∖ {i}}$ . We then average over the set of all test instances $T$ to compute the mean percentile rank (MPR):

MPR = \frac{1}{| T |} \sum Y \in T {PR}_{i, Y ∖ {i}} .

Appendix B Hyperparameters for experiments in tab:predictive-qual

Preventing numerical instabilities: The first term on the right side of eq:nonsymm-full-log-likelihood will be singular whenever $| Y_{i} | > K$ , where $Y_{i}$ is an observed subset. Therefore, to address this in practice we set $K$ to the size of the largest subset observed in the data, $K^{'}$ , as in [13]. However, this does not entirely address the issue, as the first term on the right side of eq:nonsymm-full-log-likelihood may still be singular even when $| Y_{i} | \leq K$ . In this case though, we know that we are not at a maximum, since the value of the objective function is $- \infty$ . Numerically, to prevent such singularities, in our implementation we add a small $ϵ I$ correction to each $L_{Y_{i}}$ when optimizing eq:nonsymm-full-log-likelihood (we set $ϵ = 10^{- 5}$ in our experiments).

We perform a grid search using a held-out validation set to select the best performing hyperparameters for each model and dataset. The hyperparameter settings used for each model and dataset are described below.

Symmetric low-rank DPP [12]. For this model, we use $K$ for the number of item feature dimensions for the symmetric component $V$ , and $α$ for the regularization hyperparameter for $V$ . We use the following hyperparameter settings:

Both Amazon datasets: $K = 30, α = 0$ .
UK Retail dataset: $K = 100, α = 1$ .
Instacart dataset: $K = 100, α = 0.01$ .

Baseline NDPP [11]. For this model, to ensure consistency with the notation used in [11], we use $D$ to denote the number of item feature dimensions for the symmetric component $V$ , and $D^{'}$ to denote the number of item feature dimensions for the nonsymmetric components, $B$ and $C$ . As described in [11], $α$ is the regularization hyperparameter for the $V$ , while $β$ and $γ$ are the regularization hyperparameters for $B$ and $C$ , respectively. We use the following hyperparameter settings:

Both Amazon datasets: $D = 30, α = 0$ .
Amazon apparel dataset: $D^{'} = 30$ .
Amazon three-category dataset: $D^{'} = 100$ .
UK Retail dataset: $D = 100, D^{'} = 20, α = 1$ .
All datasets: $β = γ = 0$ .

Scalable NDPP. As described in sec:model, we use $K$ to denote the number of item feature dimensions for the symmetric component $V$ and the dimensionality of the nonsymmetric component $C$ . $α$ is the regularization hyperparameter. We use the following hyperparameter settings:

Amazon apparel dataset: $K = 30, α = 0$ .
Amazon three-category dataset: $K = 80, α = 1$ .
UK dataset: $K = 100, α = 0.001$ .
Instacart dataset: $K = 100, α = 0.01$ .

For all of the above model configurations we use a batch size of 200 during training, except for the scalable NDPPs trained on the Amazon apparel and Instacart datasets, where a batch size of 400 is used.

Appendix C Training time

In fig:training-time, we report the wall-clock training time of the decomposition of [11] (NDPP) and our scalable NDPP for the Amazon: 3-category (fig:training-time:amazonthree) and UK Retail (fig:training-time:ukretail) datasets. For the Amazon: 3-category dataset, both approaches show comparable results, with the scalable NDPP converging $1.07$ times faster than NDPP. But for the UK Retail dataset, which has a much larger ground set, our scalable NDPP achieves convergence about $8.31$ times faster. We do not have a timing comparison for the Instacart dataset because the model with the decomposition of [11] cannot be trained on this dataset due to prohibitive memory and computational costs.

Appendix D Benchmark algorithms for MAP inference

We test following approximate algorithms for MAP inference:

Greedy local search.

This algorithm starts from the output of greedy, $Y^{G}$ , and replaces $i \in Y^{G}$ with $j \notin Y^{G}$ that gives the maximum improvement of the determinant, if such $i, j$ exist. Kathuria and Deshpande [22] showed that running the search for such a swap $O (k^{2} log (k / ε))$ times with an accuracy parameter $ε$ gives a tight approximation guarantee for MAP inference for symmetric DPPs. We set the number of swaps to $⌊ k^{2} log (10 k) ⌋$ for $ε = 0.1$ and use greedy local search as a baseline, since it is strictly an improvement on the greedy solution. However, greedy local search requires $O (M K k^{5} log (k / ε))$ operations, and thus it is the slowest among all of our baseline algorithms.

Stochastic greedy.

This algorithm computes marginal gains of a few items chosen uniformly at random and selects the best among them. [36] proved that $O (M / k log (1 / ε))$ samples are enough to guarantee an $(1 - 1 / e - ε)$ approximation ratio for submodular functions (i.e., symmetric DPPs). We choose $ε = 0.1$ and set the number of samples to $⌊ (M / k) log (10) ⌋$ . Under this setting, the time complexity of stochastic greedy is $O (M K k^{2} log (1 / ε))$ , which is better than the naïve exact greedy algorithm. However, we note that it is worse than that of our efficient greedy implement (alg:greedylowrank). This is because the stochastic greedy uses different random samples for every iteration and this does not take advantage of the amortized computations in lmm:bilinear-inverse. In our experiments, we simply modify line 10 in alg:greedylowrank for stochastic greedy ( $argmax$ is operated on a random subset of marginal gains), hence it can run in $O (M K k + M / k log (1 / ε))$ time. In practice, we observe that stochastic greedy is slightly slower than exact greedy due to the additional costs of random sampling process.

MCMC sampling.

We also compare inference algorithms with sampling from a DPP with fixed size (known as a $k$ -DPP). Exact sampling [27] requires eigendecomposition of $L$ , which is infeasible for a large $M$

. To resolve this, Markov Chain Monte-Carlo (MCMC) based sampling is preferred. In particular, we consider the Metropolis-Hastings algorithm, which begins with a random subset

Y

with size

k

, and picks

i \in Y

and

j \notin Y

uniformly at random. Then, it swaps them with probability

min (1, \frac{det (L_{Y \cup {j} ∖ {i}})}{det (L_{Y})})

(11)

and repeats this process for several steps. Recent work [29] shows that the MCMC sampling provides promising results for kernel approximation. We set the number of swaps to $⌊ k^{2} log (10 k) ⌋$ (the same as for greedy local search), for a runtime complexity of $O (M + K k^{4} log k)$ , which is better than the greedy algorithm.

We provide the wall-clock time of the above algorithms for real-world datasets in tab:mapinferencetime. Observe that the greedy algorithm is the fastest method for all datasets except Instacart. For Instacart, MCMC sampling is faster than other approaches, but it has much larger relative errors in terms of log-determinant (see tab:mapinference), which is not suitable for our purposes.

Algorithms

Amazon: Apparel

Amazon: 3-category

UK Retail

Instacart

Greedy local search [22]

4.28 ms

7.11 ms

36.92 ms

468.50 ms

Greedy (alg:greedylowrank)

0.19 ms

0.27 ms

1.88 ms

21.86 ms

Stochastic greedy [36]

0.26 ms

0.38 ms

2.03 ms

24.40 ms

MCMC sampling [29]

6.30 ms

6.44 ms

10.88 ms

16.86 ms

Table 4: Wall-clock time (in milliseconds) of MAP inference algorithms on NDPPs learned from real-world datasets.

Appendix E Corollary of thm:greedy

thm:greedy requires the technical condition $σ_{min} > 1$ but in practice there is no particular evidence that this condition holds. While this condition can be achieved by multiplying $L$ by a constant, this leads to a (potentially large) additive term in eq:greedyratio. Here, we provide cor:greedy which excludes the $σ_{min} > 1$ assumption from thm:greedy, and quantifies this additive term.

Corollary 1.

Consider a nonsymmetric low-rank DPP $L = V V^{⊤} + B C B^{⊤}$ , where $V, B$ are of rank $K$ , and $C \in R^{K \times K}$ . Given a cardinality budget $k$ , let $σ_{min}$ and $σ_{max}$ denote the smallest and largest singular values of $L_{Y}$ for all $Y \subseteq [[M]]$ and $| Y | \leq 2 k$ . Let $κ := σ_{max} / σ_{min}$ . Then,

log det (L_{Y^{G}})

\geq \frac{4 (1 - e^{- 1 / 4})}{2 log κ + 1} log det (L_{Y^{*}}) - (1 - \frac{4 (1 - e^{- 1 / 4})}{2 log κ + 1}) k (1 - log σ_{min})

(12)

where $Y^{G}$ is the output of alg:greedylowrank and $Y^{*}$ is the optimal solution of MAP inference in eq:mapinference.

The proof of cor:greedy is provided in sec:proof-corr-greedy. Note that instead of $log (σ_{max}) / log (σ_{min})$ , cor:greedy has a $log (σ_{max} / σ_{min})$ term in the denominator.

Appendix F Proofs

f.1 Proof of lemma:skew-sym-decomp

See 1

Proof.

First, $A = {P Σ P}^{⊤}$

for some orthogonal matrix

P \in R^{M \times M}

and

$Σ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0λ1−λ100λ2\huge 0−λ20⋱0λℓ/2\huge 0−λℓ/200⋱0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠$

(see, e.g.,[41, Proposition 2.1], which is easily extended to the case when $M$

is odd).

Let $C$ be the $ℓ \times ℓ$ submatrix of $Σ$ obtained by keeping its first $ℓ$ rows and columns and let $Q = ⎛ ⎜ ⎝ \begin{matrix} I_{ℓ} 0 \end{matrix} ⎞ ⎟ ⎠$ , where $I_{ℓ}$ is the $ℓ \times ℓ$ identity matrix. Then, $Σ = Q C Q^{⊤}$ and one can write $A = P Q C Q^{⊤} P^{⊤}$ . Setting $B = P Q$ proves the lemma. ∎

f.2 Proof of theorem:loglik-complexity

See 1

Proof.

To show that this log-likelihood can be computed in time linear in $M$ , we first show that the DPP normalization term can be computed in linear time. We briefly note here that a key component in the proof of the normalization term complexity is use of the Woodbury matrix identity, which allows us to transform an $O (M^{3})$ operation into an $O (M K^{2})$ one. A similar approach was also recently adopted in normalizing flow for deep generative models [31].

First, assume that $C$ is invertible. Otherwise, replace $C$ with $ε I_{k} + C$ for some small $ε > 0$ . Then we have:

	$det (I + L)$
	$= det (I + V V^{⊤}) det (C) det (C^{- 1} + B^{⊤} (I + V V^{⊤})^{- 1} B)$		(13)
	$= det (I_{K} + V^{⊤} V) det (C) det (C^{- 1} + B^{⊤} (I + V V^{⊤})^{- 1} B)$		(14)
	$= det (I_{K} + V^{⊤} V) det (C) det (C^{- 1} + B^{⊤} (I - V (I_{K} + V^{⊤} V)^{- 1} V^{⊤}) B)$		(15)
	$= det (I_{K} + V^{⊤} V) det (C) det (C^{- 1} + B^{⊤} B - B^{⊤} V (I_{K} + V^{⊤} V)^{- 1} V^{⊤} B) .$		(16)

eq:norm-com-proof-0 follows from application of the generalized form of the matrix determinant lemma. eq:norm-com-proof-1 follows from application of the Weinstein–Aronszajn identity. eq:norm-com-proof-2 follows from application of Woodbury’s matrix identity. If we analyze the final form in eq:norm-com-proof-3, we see that it consists of three determinants of $K \times K$ matrices. Thus, once the internal matrices are computed, the cost of the determinants is $O (K^{3})$ . The most expensive internal matrix computation is $B^{⊤} (I + V V^{⊤})^{- 1} B$ . Computing $V^{⊤} V$ takes $O (M K^{2})$ operations. Computing the inverse takes $O (K^{3})$ operations, and the remaining multiplications by $B$ require $O (M K^{2})$ operations. Therefore, computing $det (I + L)$ requires $O (M K^{2} + K^{3})$ operations overall.

Having established that the normalization term in the likelihood can be computed in $O (M K^{2} + K^{3})$ time, we proceed with characterizing the complexity of the other terms in the likelihood. The first term in eq:nonsymm-full-log-likelihood consists of determinants of size $| Y_{i} |$ . Assuming that these never exceed size $K^{'}$ , each can be computed in at most $O (K^{' 3})$ time. The regularization term is a simple sum of norms that can be computed in $O (M K)$ time. Therefore, the full regularized log-likelihood can be computed in $O (M K^{2} + K^{3} + n K^{' 3})$ time. ∎

f.3 Proof of theorem:grad-complexity

See 2

Proof.

To prove that the gradient of the log-likelihood can be computed in time linear in $M$ , we begin by showing that the gradient of the DPP’s normalization term can be computed in linear time.

From eq:norm-com-proof-1 in theorem:loglik-complexity’s proof, we have:

	$Z$	$= log det (I + L)$
		$= log det (I_{K} + V^{⊤} V) + log det (C) + log det (C^{- 1} + B^{⊤} (I + V V^{⊤})^{- 1} B)$		(17)

The gradient of $Z$ has three parts: $\nabla Z = [\nabla_{V} Z, \nabla_{B} Z, \nabla_{C} Z]$ . We derive each below, making use of the standard rule for the gradient of $log det$ : for any matrix $A$ , $\frac{\partial}{\partial A_{i j}} (log det (A)) = t r (A^{- 1} \frac{\partial A}{\partial A_{i j}}) = (A^{- 1})^{⊤}$ . If we define $X = B^{⊤} (I + V V^{⊤})^{- 1} B$ , the gradients are:

$\nabla_{V} Z$	$= \nabla_{V} log det (I_{K} + V^{⊤} V) + \nabla_{V} log det (C^{- 1} + X)$	(18)
	$= 2 V (I_{K} + V^{⊤} V)^{- 1} - (I + V V^{⊤})^{- 1} B ((C^{- 1} + X)^{- 1}$
	$+ ((C^{- 1})^{⊤} + X)^{- 1}) B^{⊤} (I + V V^{⊤})^{- 1} V$	(19)
$\nabla_{B} Z$	$= \nabla$

Scalable Learning and MAP Inference for Nonsymmetric Determinantal Point Processes

Authors

Learning Nonsymmetric Determinantal Point Processes

Wasserstein Learning of Determinantal Point Processes

Learning Determinantal Point Processes by Sampling Inferred Negatives

A Scalable CUR Matrix Decomposition Algorithm: Lower Time Complexity and Tighter Bound

Linear-time Learning on Distributions with Approximate Kernel Embeddings

The Bayesian Low-Rank Determinantal Point Process Mixture Model

Deep Determinantal Point Processes

This week in AI

1 Introduction

2 Background

3 New kernel decomposition and scalable learning

Lemma 1.

Theorem 1.

Theorem 2.

4 MAP inference

Lemma 2.

4.1 Approximation guarantee for greedy NDPP MAP inference

Theorem 3.

4.2 Greedy conditioning for next-item prediction

5 Experiments

5.1 Datasets

5.2 Experimental setup and metrics

5.3 Predictive performance results for learning

5.4 Performance results for MAP inference

5.5 Performance guarantee for greedy MAP inference

6 Conclusion

Broader Impact

References

Appendix A Mean Percentile Rank

Appendix B Hyperparameters for experiments in tab:predictive-qual

Appendix C Training time

Appendix D Benchmark algorithms for MAP inference

Greedy local search.

Stochastic greedy.

MCMC sampling.

Appendix E Corollary of thm:greedy

Corollary 1.

Appendix F Proofs

f.1 Proof of lemma:skew-sym-decomp

Proof.

f.2 Proof of theorem:loglik-complexity

Proof.

f.3 Proof of theorem:grad-complexity

Proof.