Year: 2021

Abstract:
A determinantal point process (DPP) on a collection of M items is a model, parameterized by a symmetric kernel matrix, that assigns a probability to every subset of those items. Recent work shows that removing the kernel symmetry constraint, yielding nonsymmetric DPPs (NDPPs), can lead to significant predictive performance gains for machine learning applications. However, existing work leaves open the question of scalable NDPP sampling. There is only one known DPP sampling algorithm, based on Cholesky decomposition, that can directly apply to NDPPs as well. Unfortunately, its runtime is cubic in M , and thus does not scale to large item collections. In this work, we first note that this algorithm can be transformed into a linear-time one for kernels with low-rank structure. Furthermore, we develop a scalable sublinear-time rejection sampling algorithm by constructing a novel proposal distribution. Additionally, we show that imposing certain structural constraints on the NDPP kernel enables us to bound the rejection rate in a way that depends only on the kernel rank. In our experiments we compare the speed of all of these samplers for a variety of real-world tasks.

heading: 1 INTRODUCTION
text: A determinantal point process (DPP) on M items is a model, parameterized by a symmetric kernel matrix, that assigns a probability to every subset of those items. DPPs have been applied to a wide range of machine learning tasks, including stochastic gradient descent (SGD) (Zhang et al., 2017), reinforcement learning (Osogami & Raymond, 2019; Yang et al., 2020), text summarization (Dupuy & Bach, 2018), coresets (Tremblay et al., 2019), and more. However, a symmetric kernel can only capture negative correlations between items. Recent works (Brunel, 2018; Gartrell et al., 2019) have shown that using a nonsymmetric DPP (NDPP) allows modeling of positive correlations as well, which can lead to significant predictive performance gains. Gartrell et al. (2021) provides scalable NDPP kernel learning and MAP inference algorithms, but leaves open the question of scalable sampling. The only known sampling algorithm for NDPPs is the Cholesky-based approach described in Poulson (2019), which has a runtime of O(M3) and thus does not scale to large item collections.
There is a rich body of work on efficient sampling algorithms for (symmetric) DPPs, including recent works such as Derezinski et al. (2019); Poulson (2019); Calandriello et al. (2020). Key distinctions between existing sampling algorithms include whether they are for exact or approximate sampling, whether they assume the DPP kernel has some low-rank K M , and whether they sample from the space of all 2M subsets or from the restricted space of size-k subsets, so-called k-DPPs. In this work we focus on exact sampling for low-rank kernels, and provide scalable algorithms for NDPPs. Our contributions are as follows, with runtime and memory details summarized in Table 1:
• Linear-time sampling (Section 3): We show how to transform the O(M3) Choleskydecomposition-based sampler from Poulson (2019) into an O(MK2) sampler for rank-K kernels. • Sublinear-time sampling (Section 4): Using rejection sampling, we show how to leverage existing sublinear-time samplers for symmetric DPPs to implement a sublinear-time sampler for a subclass of NDPPs that we call orthogonal NDPPs (ONDPPs).
• Learning with orthogonality constraints (Section 5): We show that the scalable NDPP kernel learning of Gartrell et al. (2021) can be slightly modified to impose an orthogonality constraint, yielding the ONDPP subclass. The constraint allows us to control the rejection sampling algorithm’s rejection rate, ensuring its scalability. Experiments suggest that the predictive performance of the kernels is not degraded by this change.
For a common large-scale setting where M is 1 million, our sublinear-time sampler results in runtime that is hundreds of times faster than the linear-time sampler. In the same setting, our linear-time sampler provides runtime that is millions of times faster than the only previously known NDPP sampling algorithm, which has cubic time complexity and is thus impractical in this scenario.

heading: 2 BACKGROUND
text: Notation. We use [M ] := {1, . . . ,M} to denote the set of items 1 through M . We use IK to denote the K-by-K identity matrix, and often write I := IM when the dimensionality should be clear from context. Given L ∈ RM×M , we use Li,j to denote the entry in the i-th row and j-th column, and LA,B ∈ R|A|×|B| for the submatrix formed by taking rows A and columns B. We also slightly abuse notation to denote principal submatrices with a single subscript, LA := LA,A.
Kernels. As discussed earlier, both (symmetric) DPPs and NDPPs define a probability distribution over all 2M subsets of a ground set [M ]. The distribution is parameterized by a kernel matrix L ∈ RM×M and the probability of a subset Y ⊆ [M ] is defined to be Pr(Y ) ∝ det(LY ). For this to define a valid distribution, it must be the case that det(LY ) ≥ 0 for all Y . For symmetric DPPs, the non-negativity requirement is identical to a requirement that L be positive semi-definite (PSD). For nonsymmetric DPPs, there is no such simple correspondence, but prior work such as Gartrell et al. (2019; 2021) has focused on PSD matrices for simplicity.
Normalizing and marginalizing. The normalizer of a DPP or NDPP distribution can also be written as a single determinant: ∑ Y⊆[M ] det(LY ) = det(L+ I) (Kulesza & Taskar, 2012, Theorem 2.1). Additionally, the marginal probability of a subset can be written as a determinant: Pr(A ⊆ Y ) = det(KA), for K := I − (L+ I)−1 (Kulesza & Taskar, 2012, Theorem 2.2)*, where K is typically called the marginal kernel.
Intuition. The diagonal element Ki,i is the probability that item i is included in a set sampled from the model. The 2-by-2 determinant det(K{i,j}) = Ki,iKj,j −Ki,jKj,j is the probability that both i and j are included in the sample. A symmetric DPP has a symmetric marginal kernel, meaning Ki,j = Kj,i, and hence Ki,iKj,j −Ki,jKj,i ≤ Ki,iKj,j . This implies that the probability of including both i and j in the sampled set cannot be greater than the product of their individual inclusion probabilities. Hence, symmetric DPPs can only encode negative correlations. In contrast, NDPPs can have Ki,j and Kj,i with differing signs, allowing them to also capture positive correlations.

heading: 2.1 RELATED WORK
text: Learning. Gartrell et al. (2021) proposes a low-rank kernel decomposition for NDPPs that admits linear-time learning. The decomposition takes the form L := V V > + B(D − D>)B> for V ,B ∈ RM×K , and D ∈ RK×K . The V V > component is a rank-K symmetric matrix, which can model negative correlations between items. The B(D −D>)B> component is a rank-K skewsymmetric matrix, which can model positive correlations between items. For compactness of notation, we will write L = ZXZ>, where Z = [ V B ] ∈ RM×2K , and X = [ IK 0
0 D−D>
] ∈ R2K×2K .
The marginal kernel in this case also has a rank-2K decomposition, as can be shown via application of the Woodbury matrix identity:
K := I − (I +L)−1 = ZX ( I2K +Z >ZX )−1 Z>. (1)
*The proofs in Kulesza & Taskar (2012) typically assume a symmetric kernel, but this particular one does not rely on the symmetry.
Algorithm 1 Cholesky-based NDPP sampling (Poulson, 2019, Algorithm 1) 1: procedure SAMPLECHOLESKY(K) . marginal kernel factorization Z,W 2: Y ← ∅ Q←W 3: for i = 1 to M do 4: pi ←Ki,i pi ← z>i Qzi 5: u← uniform(0, 1) 6: if u ≤ pi then Y ← Y ∪ {i} 7: else pi ← pi − 1 8: KA ←KA − KA,iKi,Api for A := {i+ 1, . . . ,M} Q← Q− Qziz > i Q pi
9: return Y
Note that the matrix to be inverted can be computed from Z and X in O(MK2) time, and the inverse itself takes O(K3) time. Thus, K can be computed from L in time O(MK2). We will develop sampling algorithms for this decomposition, as well as an orthogonality-constrained version of it. We use W := X ( I2K +Z >ZX )−1 in what follows so that we can compactly write K = ZWZ>.
Sampling. While there are a number of exact sampling algorithms for DPPs with symmetric kernels, the only published algorithm that clearly can directly apply to NDPPs is from Poulson (2019) (see Theorem 2 therein). This algorithm begins with an empty set Y = ∅ and iterates through the M items, deciding for each whether or not to include it in Y based on all of the previous inclusion/exclusion decisions. Poulson (2019) shows, via the Cholesky decomposition, that the necessary conditional probabilities can be computed as follows:
Pr (j ∈ Y | i ∈ Y ) = Pr({i, j} ⊆ Y ) Pr(i ∈ Y ) = Kj,j − (Kj,iKi,j) /Ki,i, (2)
Pr (j ∈ Y | i /∈ Y ) = Pr(j ∈ Y )− Pr({i, j} ⊆ Y ) Pr(i /∈ Y ) = Kj,j − (Kj,iKi,j) / (Ki,i − 1) . (3)
Algorithm 1 (left-hand side) gives pseudocode for this Cholesky-based sampling algorithm†.
There has also been some recent work on approximate sampling for fixed-size k-NDPPs: Alimohammadi et al. (2021) provide a Markov chain Monte Carlo (MCMC) algorithm and prove that the overall runtime to approximate ε-close total variation distance is bounded by O(M2k3 log(1/(εPr(Y0))), where Pr(Y0) is probability of an initial state Y0. Improving this runtime is an interesting avenue for future work, but for this paper we focus on exact sampling.

heading: 3 LINEAR-TIME CHOLESKY-BASED SAMPLING
text: In this section, we show that the O(M3) runtime of the Cholesky-based sampler from Poulson (2019) can be significantly improved when using the low-rank kernel decomposition of Gartrell et al. (2021). First, note that Line 8 of Algorithm 1, where all marginal probabilities are updated via an (M − i)-by-(M − i) matrix subtraction, is the most costly part of the algorithm, making overall time and memory complexities O(M3) and O(M2), respectively. However, when the DPP kernel is given by a low-rank decomposition, we observe that marginal probabilities can be updated by matrix-vector multiplications of dimension 2K, regardless of M . In more detail, suppose we have the marginal kernel K = ZWZ> as in Eq. (1) and let zj be the j-th row vector in Z. Then, for i 6= j:
Pr (j ∈ Y | i ∈ Y ) = Kj,j − (Kj,iKi,j)/Ki,i = z>j ( W − (Wzi)(z > i W )
z>i Wzi
) zj , (4)
Pr (j ∈ Y | i /∈ Y ) = z>j ( W − (Wzi)(z > i W )
z>i Wzi − 1
) zj . (5)
The conditional probabilities in Eqs. (4) and (5) are of bilinear form, and the zj do not change during sampling. Hence, it is enough to update the 2K-by-2K inner matrix at each iteration, and obtain the marginal probability by multiplying this matrix by zi. The details are shown on the right-hand side of Algorithm 1. The overall time and memory complexities are O(MK2) and O(MK), respectively.
†Cholesky decomposition is defined only for a symmetric positive definite matrix. However, we use the term “Cholesky” from Poulson (2019) to maintain consistency with this work, although Algorithm 1 is valid for nonsymmetric matrices.

heading: 4 SUBLINEAR-TIME REJECTION SAMPLING
text: Although the Cholesky-based sampler runs in time linear in M , even this is too expensive for the large M that are often encountered in real-world datasets. To improve runtime, we consider rejection sampling (Von Neumann, 1963). Let p be the target distribution that we aim to sample, and let q be any distribution whose support corresponds to that of p; we call q the proposal distribution. Assume that there is a universal constant U such that p(x) ≤ Uq(x) for all x. In this setting, rejection sampling draws a sample x from q and accepts it with probability p(x)/(Uq(x)), repeating until an acceptance occurs. The distribution of the resulting samples is p. It is important to choose a good proposal distribution q so that sampling is efficient and the number of rejections is small.
Algorithm 2 Rejection NDPP sampling (Tree-based sampling) 1: procedure PREPROCESS(V ,B,D) 2: {(σj ,y2j−1,y2j)}K/2j=1 ← YOULADECOMPOSE(B,D)‡ 3: X̂ ← diag ( IK , σ1, σ1, . . . , σK/2, σK/2
) 4: Z ← [V ,y1, . . . ,yK ] {(λi, zi)}2Ki=1 ← EIGENDECOMPOSE(ZX̂1/2) T ← CONSTRUCTTREE(M, [z1, . . . ,z2K ]>) 5: return Z, X̂ return T , {(λi, zi)}2Ki=1 6: procedure SAMPLEREJECT(V ,B,D,Z, X̂) . tree T , eigen pair {(λi, zi)}2Ki=1 of ZX̂Z 7: while true do 8: Y ← SAMPLEDPP(ZX̂Z>) Y ← SAMPLEDPP(T , {(λi, zi)}2Ki=1) 9: u← uniform(0, 1)
10: p← det([V V >+B(D−D>)B>]Y )
det([ZX̂Z>]Y )
11: if u ≤ p then break 12: return Y

heading: 4.1 PROPOSAL DPP CONSTRUCTION
text: Our first goal is to find a proposal DPP with symmetric kernel L̂ that can upper-bound all probabilities of samples from the NDPP with kernel L within a constant factor. To this end, we expand the determinant of a principal submatrix, det(LY ), using the spectral decomposition of the NDPP kernel. Such a decomposition essentially amounts to combining the eigendecomposition of the symmetric part of L with the Youla decomposition (Youla, 1961) of the skew-symmetric part.
Specifically, suppose {(σj ,y2j−1,y2j)}K/2j=1 is the Youla decomposition of B(D −D>)B> (see Appendix D for more details), that is,
B(D −D>)B> = K/2∑ j=1 σj ( y2j−1y > 2j − y2jy>2j−1 ) . (6)
Then we can simply write L = ZXZ>, for Z := [V ,y1, . . . ,yK ] ∈ RM×2K , and
X := diag ( IK , [ 0 σ1 −σ1 0 ] , . . . , [ 0 σK/2 −σK/2 0 ]) . (7)
Now, consider defining a related but symmetric PSD kernel L̂ := ZX̂Z> with X̂ := diag ( IK , σ1, σ1, . . . , σK/2, σK/2 ) . All determinants of the principal submatrices of L̂ = ZX̂Z> upper-bound those of L, as stated below.
Theorem 1. For every subset Y ⊆ [M ], it holds that det(LY ) ≤ det(L̂Y ). Moreover, equality holds when the size of Y is equal to the rank of L.
‡Pseudo-code of YOULADECOMPOSE is provided in Algorithm 4. See Appendix D.
Algorithm 3 Tree-based DPP sampling (Gillenwater et al., 2019) 1: procedure BRANCH(A,Z) 2: if A = {j} then 3: T .A← {j}, T .Σ← Z>j,:Zj,: 4: return T 5: A`, Ar ← Split A in half 6: T .left← BRANCH(A`,Z) 7: T .right← BRANCH(Ar,Z) 8: T .Σ← T .left.Σ + T .right.Σ 9: return T 10: procedure CONSTRUCTTREE(M , Z) 11: return BRANCH([M ], Z)
12: procedure SAMPLEDPP(T ,Z, {λi}Ki=1) 13: E ← ∅, Y ← ∅, QY ← 0 14: for i = 1, . . . ,K do 15: E ← E ∪ {i} w.p. λi/(λi + 1) 16: for k = 1, . . . , |E| do 17: j ← SAMPLEITEM(T ,QY , E) 18: Y ← Y ∪ {j} 19: QY← I|E|−Z>Y,E ( ZY,EZ > Y,E )−1 ZY,E
20: return Y
21: procedure SAMPLEITEM(T ,QY , E) 22: if T is a leaf then return T .A 23: p` ← 〈 T .left.ΣE ,QY 〉 24: pr ← 〈 T .right.ΣE ,QY
〉 25: u← uniform(0, 1) 26: if u ≤ p`p`+pr then 27: return SAMPLEITEM(T .left,QY , E) 28: else 29: return SAMPLEITEM(T .right,QY , E)
Proof sketch: From the Cauchy-Binet formula, the determinants of LY and L̂Y for all Y ⊆ [M ], |Y | ≤ 2K can be represented as
det(LY ) = ∑
I⊆[K],|I|=|Y | ∑ J⊆[K],|J|=|Y | det(XI,J) det(ZY,I) det(ZY,J), (8)
det(L̂Y ) = ∑
I⊆[2K],|I|=|Y |
det(X̂I) det(ZY,I) 2. (9)
Many of the terms in Eq. (8) are actually zero due to the block-diagonal structure of X . For example, note that if 1 ∈ I but 1 /∈ J , then there is an all-zeros row in XI,J , making det(XI,J) = 0. We show that each XI,J with nonzero determinant is a block-diagonal matrix with diagonal entries among
±σj , or [ 0 σj −σj 0 ] . With this observation, we can prove that det(XI,J) is upper-bounded by det(X̂I)
or det(X̂J). Then, through application of the rearrangement inequality, we can upper-bound the sum of the det(XI,J) det(ZY,I) det(ZY,J) in Eq. (8) with a sum over det(X̂I) det(ZY,I)2. Finally, we show that the number of non-zero terms in Eq. (8) is identical to the number of non-zero terms in Eq. (9). Combining these gives us the desired inequality det(LY ) ≤ det(L̂Y ). The full proof of Theorem 1 is in Appendix E.1.
Now, recall that the normalizer of a DPP (or NDPP) with kernel L is det(L + I). The ratio of probability of the NDPP with kernel L to that of a DPP with kernel L̂ is thus:
PrL(Y ) PrL̂(Y ) = det(LY )/det(L+ I) det(L̂Y )/det(L̂+ I) ≤ det(L̂+ I) det(L+ I) ,
where the inequality follows from Theorem 1. This gives us the necessary universal constant U upper-bounding the ratio of the target distribution to the proposal distribution. Hence, given a sample Y drawn from the DPP with kernel L̂, we can use acceptance probability PrL(Y )/(U PrL̂(Y )) = det(LY )/ det(L̂Y ). Pseudo-codes for proposal construction and rejection sampling are given in Algorithm 2. Note that to derive L̂ from L it suffices to run the Youla decomposition of B(D − D>)B>, because the difference is only in the skew-symmetric part. This decomposition can run in O(MK2) time; more details are provided in Appendix D. Since L̂ is a symmetric PSD matrix, we can apply existing fast DPP sampling algorithms to sample from it. In particular, in the next section we combine a fast tree-based method with rejection sampling.

heading: 4.2 SUBLINEAR-TIME TREE-BASED SAMPLING
text: There are several DPP sampling algorithms that run in sublinear time, such as tree-based (Gillenwater et al., 2019) and intermediate (Derezinski et al., 2019) sampling algorithms. Here, we consider applying the former, a tree-based approach, to sample from the proposal distribution defined by L̂.
We give some details of the sampling procedure, as in the course of applying it we discovered an optimization that slightly improves on the runtime of prior work. Formally, let {(λi, zi)}2Ki=1 be the eigendecomposition of L̂ and Z := [z1, . . . ,z2K ] ∈ RM×2K . As shown in Kulesza & Taskar (2012, Lemma 2.6), for every Y ⊆ [M ], |Y | ≤ 2K, the probability of Y under DPP with L̂ can be written:
PrL̂(Y ) = det(L̂Y )
det(L̂+ I) = ∑ E⊆[2K],|E|=|Y | det(ZY,EZ > Y,E) ∏ i∈E λi λi + 1 ∏ i/∈E 1 λi + 1 . (10)
A matrix of the form Z:,EZ>:,E can be a valid marginal kernel for a special type of DPP, called an elementary DPP. Hence, Eq. (10) can be thought of as DPP probabilities expressed as a mixture of elementary DPPs. Based on this mixture view, DPP sampling can be done in two steps: (1) choose an elementary DPP according to its mixture weight, and then (2) sample a subset from the selected elementary DPP. Step (1) can be performed by 2K independent random coin tossings, while step (2) involves computational overhead. The key idea of tree-based sampling is that step (2) can be accelerated by traversing a binary tree structure, which can be done in time logarithmic in M .
More specifically, given the marginal kernel K = Z:,EZ>:,E , where E is obtained from step (1), we start from the empty set Y = ∅ and repeatedly add an item j to Y with probability: Pr(j ∈ S | Y ⊆ S) = Kj,j −Kj,Y (KY )−1KY,j = Zj,EQYZ>j,E = 〈 QY , (Z>j,:Zj,:)E 〉 , (11)
where S is some final selected subset, and QY := I|E| − Z>Y,E ( ZY,EZ > Y,E )−1 ZY,E . Consider a binary tree whose root includes a ground set [M ]. Every non-leaf node contains a subsetA ⊆ [M ] and stores a 2K-by-2K matrix ∑ j∈AZ > j,:Zj,:. A partitionA` andAr, such thatA`∪Ar = A,A`∩Ar = ∅, are passed to its left and right subtree, respectively. The resulting tree has M leaves and each has exactly a single item. Then, one can sample a single item by recursively moving down to the left node with probability:
p` = 〈QY ,∑j∈A`(Z>j,:Zj,:)E〉 〈QY ,∑j∈A(Zj,:Z>j,:)E〉 , (12)
or to the right node with probability 1− p`, until reaching a leaf node. An item in the leaf node is chosen with probability according to Eq. (11). Since every subset in the support of an elementary DPP with a rank-k kernel has exactly k items, this process is repeated for |E| iterations. Full descriptions of tree construction and sampling are provided in Algorithm 3. The proposed tree-based rejection sampling for an NDPP is outlined on the right-side of Algorithm 2. The one-time pre-processing step of constructing the tree (CONSTRUCTTREE) requires O(MK2) time. After pre-processing, the procedure SAMPLEDPP involves |E| traversals of a tree of depth O(logM), where in each node a O(|E|2) operation is required. The overall runtime is summarized in Proposition 1 and the proof can be found in Appendix E.2. Proposition 1. The tree-based sampling procedure SAMPLEDPP in Algorithm 3 runs in time O(K + k3 logM + k4), where k is the size of the sampled set§.

heading: 4.3 AVERAGE NUMBER OF REJECTIONS
text: We now return to rejection sampling and focus on the expected number of rejections. The number of rejections of Algorithm 2 is known to be a geometric random variable with mean equal to the constant U used to upper-bound the ratio of the target distribution to the proposal distribution: det(L̂+ I)/ det(L+ I). If all columns in V and B are orthogonal, which we denote V ⊥ B, then the expected number of rejections depends only on the eigenvalues of the skew-symmetric part of the NDPP kernel.
§Computing p` via Eq. (12) improves on Gillenwater et al. (2019)’s O(k4 logM) runtime for this step.
Theorem 2. Given an NDPP kernel L = V V > + B(D −D>)B> for V ,B ∈ RM×K ,D ∈ RK×K , consider the proposal kernel L̂ as proposed in Section 4.1. Let {σj}K/2j=1 be the positive eigenvalues obtained from the Youla decomposition of B(D−D>)B>. If V ⊥ B, then det(L̂+I)det(L+I) =∏K/2 j=1 ( 1 +
2σj σ2j+1 ) ≤ (1 + ω)K/2, where ω = 2K ∑K/2 j=1 2σj σ2j+1 ∈ (0, 1].
Proof sketch: Orthogonality between V and B allows det(L+ I) to be expressed just in terms of the eigenvalues of V V > and B(D −D>)B>. Since both L and L̂ share the symmetric part V V >, the ratio of determinants only depends on the skew-symmetric part. A more formal proof appears in Appendix E.3.
Assuming we have a kernel where V ⊥ B, we can combine Theorem 2 with the tree-based rejection sampling algorithm (right-side in Algorithm 2) to sample in timeO((K+k3 logM+k4)(1+ω)K/2). Hence, we have a sampling algorithm that is sublinear inM , and can be much faster than the Choleskybased algorithm when (1 + ω)K/2 M . In the next section, we introduce a learning scheme with the V ⊥ B constraint, as well as regularization to ensure that ω is small.

heading: 5 LEARNING WITH ORTHOGONALITY CONSTRAINTS
text: We aim to learn a NDPP that provides both good predictive performance and a low rejection rate. We parameterize our NDPP kernel matrix L = V V > +B(D −D>)B> by
D = diag ([ 0 σ1 0 0 ] , . . . , [ 0 σK/2 0 0 ]) (13)
for σj ≥ 0, B>B = I , and, motivated by Theorem 2, require V >B = 0¶. We call such orthogonality-constrained NDPPs “ONDPPs”. Notice that if V ⊥ B, then L has the full rank of 2K, since the intersection of the column spaces spanned by V and by B is empty, and thus the full rank available for modeling can be used. Thus, this constraint can also be thought of as simply ensuring that ONDPPs use the full rank available to them.
Given example subsets {Y1, . . . , Yn} as training data, learning is done by minimizing the regularized negative log-likelihood:
min V ,B,{σj}K/2j=1
− 1 n n∑ i=1 log ( det(LYi) det(L+ I) ) + α M∑ i=1 ‖vi‖22 µi + β M∑ i=1 ‖bi‖22 µi + γ K/2∑ j=1 log ( 1 + 2σj σ2j + 1 ) (14)
where α, β, γ > 0 are hyperparameters, µi is the frequency of item i in the training data, and vi and bi represent the rows of V and B, respectively. This objective is very similar to that of Gartrell et al. (2021), except for the orthogonality constraint and the final regularization term. Note that this regularization term corresponds exactly to the logarithm of the average rejection rate, and therefore should help to control the number of rejections.

heading: 6 EXPERIMENTS
text: We first show that the orthogonality constraint from Section 5 does not degrade the predictive performance of learned kernels. We then compare the speed of our proposed sampling algorithms.

heading: 6.1 PREDICTIVE PERFORMANCE RESULTS FOR NDPP LEARNING
text: We benchmark various DPP models, including symmetric (Gartrell et al., 2017), nonsymmetric for scalable learning (Gartrell et al., 2021), as well as our ONDPP kernels with and without rejection rate
¶Technical details: To learn NDPP models with the constraint V >B = 0, we project V according to: V ← V −B(B>B)−1(B>V ). For the B>B = I constraint, we apply QR decomposition on B. Note that both operations require O(MK2) time. (Constrained learning code is provided in the supplement. We use Pytorch’s linalg.solve to avoid the expense of explicitly computing the (B>B)−1 inverse.) Hence, our learning time complexity is identical to that of Gartrell et al. (2021).
regularization. We use the scalable NDPP models (Gartrell et al., 2021) as a baseline||. The kernel components of each model are learned using five real-world recommendation datasets, which have ground set sizes that range from 3,941 to 1,059,437 items (see Appendix A for more details).
Our experimental setup and metrics mirror those of Gartrell et al. (2021). We report the mean percentile rank (MPR) metric for a next-item prediction task, the AUC metric for subset discrimination, and the log-likelihood of the test set; see Appendix B for more details on the experiments and metrics. For all metrics, higher numbers are better. For NDPP models, we additionally report the average rejection rates when they apply to rejection sampling.
In Table 2, we observe that the predictive performance of our ONDPP models generally match or sometimes exceed the baseline. This is likely because the orthogonality constraint enables more effective use of the full rank-2K feature space. Moreover, imposing the regularization on rejection rate, as shown in Eq. (14), often leads to dramatically smaller rejection rates, while the impact on predictive performance is generally marginal. These results justify the ONDPP and regularization for fast sampling. Finally, we observe that the learning time of our ONDPP models is typically a bit longer than that of the NDPP models, but still quite reasonable (e.g., the time per iteration for the NDPP takes 27 seconds for the Book dataset, while our ONDPP takes 49.7 seconds).
Fig. 1 shows how the regularizer γ affects the test log-likelihood and the average number of rejections. We see that γ degrades predictive performance and reduces the rejection rate when set above a certain threshold; this behavior is seen for many datasets. However, for the Recipe dataset we observed that the test log-likelihood is not very sensitive to γ, likely because all models in our experiments achieve very high performance on this dataset. In general, we observe that γ can be set to a value that results in a small rejection rate, while having minimal impact on predictive performance.

heading: 6.2 SAMPLING TIME COMPARISON
text: We benchmark the Cholesky-based sampling algorithm (Algorithm 1) and tree-based rejection sampling algorithm (Algorithm 2) on ONDPPs with both synthetic and real-world data.
Synthetic datasets. We generate non-uniform random features for V ,B as done by (Han & Gillenwater, 2020). In particular, we first sample x1, . . . ,x100 from N (0, I2K/(2K)), and integers t1, . . . , t100 from Poisson distribution with mean 5, rescaling the integers such that ∑ i ti =M . Next, we draw ti random vectors from N (xi, I2K), and assign the first K-dimensional vectors as the row ||We use the code from https://github.com/cgartrel/scalable-nonsymmetric-DPPs for the NDPP baseline, which is made available under the MIT license. To simplify learning and MAP inference, Gartrell et al. (2021) set B = V in their experiments. However, since we have the V ⊥ B constraint in our ONDPP approach, we cannot set B = V . Hence, for a fair comparison, we do not set B = V for the NDPP baseline in our experiments, and thus the results in Table 2 differ slightly from those published in Gartrell et al. (2021).
10−8 10−6 10−4 10−2 100
regularizer γ
101
104
107
1010
av er
ag e
# of
re je
ct io
ns
(a)
10−6 10−4 10−2 100
regularizer γ
−106
−104
−102
−100
−98
te st
lo g-
lik el
ih oo
d
(b)
Figure 1: Average number of rejections and test log-likelihood with different values of the regularizer γ for ONDPPs trained on the UK Retail dataset. Shaded regions are 95% confidence intervals of 10 independent trials.
212 214 216 218 220
ground set size M
101
102
103
sa m
pl in
g tim
e (s
ec ) Cholesky-based Rejection
(a)
212 214 216 218 220
ground set size M
10−1
101
103
pr ep
ro ce
ss in
g tim
e (s
ec )
Tree construction Spectral decomposition
(b)
Figure 2: Wall-clock time (sec) for synthetic data for (a) NDPP sampling algorithms and (b) preprocessing steps for the rejection sampling. Shaded regions are 95% confidence intervals from 100 independent trials.
vectors of V and the latter vectors as those of B. Each entry of D is sampled from N (0, 1). We choose K = 100 and vary M from 212 to 220.
Fig. 2(a) illustrates the runtimes of Algorithms 1 and 2. We verify that the rejection sampling time tends to increase sub-linearly with the ground set size M , while the Cholesky-based sampler runs in linear time. In Fig. 2(b), the runtimes of the preprocessing steps for Algorithm 2 (i.e., spectral decomposition and tree construction) are reported. Although the rejection sampler requires these additional processes, they are one-time steps and run much faster than a single run of the Choleksy-based method for M = 220.
Real-world datasets. In Table 3, we report the runtimes and speedup of NDPP sampling algorithms for real-world datasets. All NDPP kernels are obtained using learning with orthogonality constraints, with rejection rate regularization as reported in Section 6.1. We observe that the tree-based rejection sampling runs up to 246 times faster than the Cholesky-based algorithm. For larger datasets, we expect that this gap would significantly increase. As with the synthetic experiments, we see that the tree construction pre-processing time is comparable to the time required to draw a single sample via the other methods, and thus the tree-based method is often the best choice for repeated sampling**.

heading: 7 CONCLUSION
text: In this work we developed scalable sampling methods for NDPPs. One limitation of our rejection sampler is its practical restriction to the ONDPP subclass. Other opportunities for future work include the extension of our rejection sampling approach to the generation of fixed-size samples (from k-NDPPs), the development of approximate sampling techniques, and the extension of DPP samplers along the lines of Derezinski et al. (2019); Calandriello et al. (2020) to NDPPs. Scalable sampling also opens the door to using NDPPs as building blocks in probabilistic models.
**We note that the tree can consume substantial memory, e.g., 169.5 GB for the Book dataset with K = 100. For settings where this scale of memory use is unacceptable, we suggest use of the intermediate sampling algorithm (Calandriello et al., 2020) in place of tree-based sampling. The resulting sampling algorithm may be slower, but the O(M +K) memory cost is substantially lower.

heading: A FULL DETAILS ON DATASETS
text: We perform experiments on several real-world public datasets composed of subsets:
• UK Retail: This dataset (Chen et al., 2012) contains baskets representing transactions from an online retail company that sells 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 M = 3,941 products. Baskets containing more than 100 items are in the long tail of the basket-size distribution, so omitting these is reasonable, and allows us to use a low-rank factorization of the NDPP with K = 100. • Recipe: This dataset (Majumder et al., 2019) contains recipes and food reviews from Food.com (formerly Genius Kitchen)††. Each recipe (“basket”) is composed of a collection of ingredients, resulting in 178,265 recipes and a catalog of 7,993 ingredients. • Instacart: This dataset (Instacart, 2017) contains baskets purchased by Instacart users‡‡. We omit baskets with more than 100 items, resulting in 3.2 million baskets and a catalog of 49,677 products. • Million Song: This dataset (McFee et al., 2012) contains playlists (“baskets”) of songs from Echo Nest users§§. We trim playlists with more than 100 items, leaving 968,674 playlists and a catalog of 371,410 songs. • Book: This dataset (Wan & McAuley, 2018) contains reviews from the Goodreads book review website, including a variety of attributes describing the items¶¶. For each user we build a subset (“basket”) containing the books reviewed by that user. We trim subsets with more than 100 books, resulting in 430,563 subsets and a catalog of 1,059,437 books.
As far as we are aware, these datasets do not contain personally identifiable information or offensive content. While the UK Retail dataset is publicly available, we were unable to find a license for it. Also, we were not able to determine if user consent was explicitly obtained by the organizations that constructed these datasets.

heading: B FULL DETAILS ON EXPERIMENTAL SETUP AND METRICS
text: We use 300 randomly-selected baskets as a held-out validation set, for tracking convergence during training and for tuning hyperparameters. Another 2000 random 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 & Ba, 2015) for optimization. We initialize D from the standard Gaussian distribution N (0, 1), while V and B are initialized from the uniform(0, 1) distribution.
Subset expansion task. We use greedy conditioning to do next-item prediction (Gartrell et al., 2021, Section 4.2). 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 Appendix B.1 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 to the subset.
††See https://www.kaggle.com/shuyangli94/food-com-recipes-and-user-interactions for the license for this public dataset.
‡‡This public dataset is available for non-commercial use; see https://www.instacart.com/datasets/ grocery-shopping-2017 for the license.
§§See http://millionsongdataset.com/faq/ for the license for this public dataset. ¶¶This public dataset is available for academic use only; see https://sites.google.com/eng.ucsd.edu/
ucsdbookgraph/home for the license.
B.1 MEAN PERCENTILE RANK
We begin our definition of MPR by defining percentile rank (PR). First, given a set J , let pi,J = Pr(J ∪ {i} | J). The percentile rank of an item i given a set J is defined as
PRi,J =
∑ i′ 6∈J 1(pi,J ≥ pi′,J)
|Y\J | × 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 ∈ Y and compute PRi,Y \{i}. We then average over the set of all test instances T to compute the mean percentile rank (MPR):
MPR = 1 |T | ∑ Y ∈T PRi,Y \{i}.

heading: C HYPERPARAMETERS FOR EXPERIMENTS
text: Preventing numerical instabilities: The det(LYi) in Eq. (14) will be zero whenever |Yi| > K, where Yi is an observed subset. To address this in practice we set K to the size of the largest subset observed in the data, K ′, as in Gartrell et al. (2017). However, this does not entirely fix the issue, as there is still a chance that the term will be zero even when |Yi| ≤ K. In this case though, we know that we are not at a maximum, since the value of the objective function is −∞. Numerically, to prevent such singularities, in our implementation we add a small I correction to each LYi when optimizing Eq. (14) ( = 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 (Gartrell et al., 2017). 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:
• UK Retail dataset: K = 100, α = 1. • Recipe dataset: K = 100, α = 0.01 • Instacart dataset: K = 100, α = 0.001. • Million Song dataset: K = 100, α = 0.0001. • Book dataset: K = 100, α = 0.001
Scalable NDPP (Gartrell et al., 2021). As described in Section 2.1, we useK to denote the number of item feature dimensions for the symmetric component V and the dimensionality of the nonsymmetric component D. α and β are the regularization hyperparameters. We use the following hyperparameter settings:
• UK dataset: K = 100, α = 0.01. • Recipe dataset: K = 100, α = β = 0.01. • Instacart dataset: K = 100, α = 0.001. • Million Song dataset: K = 100, α = 0.01. • Book dataset: K = 100, α = β = 0.1
ONDPP. As described in Section 5, we use K to denote the number of item feature dimensions for the symmetric component V and the dimensionality of the nonsymmetric component C. α, β, and γ are the regularization hyperparameters. We use the following hyperparameter settings:
• UK dataset: K = 100, α = β = 0.01, γ = 0.5. • Recipe dataset: K = 100, α = β = 0.01, γ = 0.1. • Instacart dataset: K = 100, α = β = 0.001, γ = 0.001. • Million Song dataset: K = 100, α = β = 0.01, γ = 0.2.
• Book dataset: K = 100, α = β = 0.01, γ = 0.1.
For all of the above model configurations and datasets, we use a batch size of 800 during training.

heading: D YOULA DECOMPOSITION: SPECTRAL DECOMPOSITION FOR SKEW-SYMMETRIC MATRIX
text: We provide some basic facts on the spectral decomposition of a skew-symmetric matrix, and introduce an efficient algorithm for this decomposition when it is given by a low-rank factorization. We write i := √ −1 and vH as the conjugate transpose of v ∈ CM , and denote Re(z) and Im(z) by the real and imaginary parts of a complex number z, respectively.
Given B ∈ RM×K and D ∈ RK×K , consider a rank-K skew-symmetric matrix B(D −D>)B>. Note that all nonzero eigenvalues of a real-valued skew-symmetric matrix are purely imaginary. Denote iσ1,−iσ1, . . . , iσK/2,−iσK/2 by its nonzero eigenvalues where each of σj is real, and a1 + ib1,a1 − ib1, . . .aK/2 + ibK/2,aK/2 − ibK/2 by the corresponding eigenvectors for aj , bj ∈ RM , which come in conjugate pairs. Then, we can write
B(D −D>)B> = K/2∑ j=1 iσj(aj + ibj)(aj + ibj) H − iσj(aj − ibj)(aj − ibj)H (15)
= K/2∑ j=1 2σj(ajb > j − bja>j ) (16)
= K/2∑ j=1 [ aj − bj aj + bj ] [ 0 σj −σj 0 ] [ a>j − b>j a>j + b > j ] . (17)
Note that a1 ± b1, . . . ,aK/2 ± bK/2 are real-valued orthonormal vectors, because a1, b1, . . . ,aK/2, bK/2 are orthogonal to each other and ‖aj ± bj‖22 = ‖aj‖ 2 2 + ‖bj‖ 2 2 = 1 for all j. The pair {(σj ,aj − bj ,aj + bj)}K/2j=1 is often called the Youla decomposition (Youla, 1961) of B(D −D>)B>. To efficiently compute the Youla decomposition of a rank-K matrix, we use the following result.
Proposition 2 (Proposition 1, Nakatsukasa (2019)). Given A,B ∈ CM×K , the nonzero eigenvalues of AB> ∈ CM×M and B>A ∈ CK×K are identical. In addition, if (λ,v) is an eigenpair of B>A with λ 6= 0, then (λ,Av/ ‖Av‖2) is an eigenpair of AB>.
From the above proposition, one can first compute (D −D>)B>B and then apply the eigendecomposition to that K-by-K matrix. Taking the imaginary part of the obtained eigenvalues gives us the σj’s, and multiplying B by the eigenvectors gives us the eigenvectors of B(D −D>)B>. In addition, this can be done in O(MK2 +K3) time; when M > K it runs much faster than the eigendecomposition of B(D −D>)B>, which requires O(M3) time. The pseudo-code of the Youla decomposition is provided in Algorithm 4.
Algorithm 4 Youla decomposition of low-rank skew-symmetric matrix 1: procedure YOULADECOMPOSITION(B,D) 2: {(ηj , zj), (ηj , zj)}K/2j=1 ← eigendecomposition of (D −D>)B>B 3: for j = 1, . . . ,K/2 do 4: σj ← Im(ηj) for j = 1, . . . ,K/2 5: y2j−1 ← B (Re(zj)− Im(zj)) 6: y2j ← B (Re(zj) + Im(zj)) 7: yj ← yj/ ‖yj‖ for j = 1, . . . ,K 8: return {(σj ,y2j−1,y2j)}K/2j=1

heading: E PROOFS
text: E.1 PROOF OF THEOREM 1
Theorem 1. For every subset Y ⊆ [M ], it holds that det(LY ) ≤ det(L̂Y ). Moreover, equality holds when the size of Y is equal to the rank of L.
Proof of Theorem 1. It is enough to fix Y ⊆ [M ] such that 1 ≤ |Y | ≤ 2K, because the rank of both L and L̂ is up to 2K. Denote k := |Y | and ( [2K] k ) := {I ⊆ [2K]; |I| = k} for k ≤ 2K. We recall the definition of L̂: given V ,B,D such that L = V V > + B(D −D>)B>, let {(ρi,vi)}Ki=1 be the eigendecomposition of V V > and {(σj ,y2j−1,y2j)}K/2j=1 be the Youla decomposition of B(D −D>)B>. Denote Z := [v1, . . . ,vK ,y1, . . . ,yK ] ∈ RM×2K and
X := diag ( ρ, . . . , ρK , [ 0 σ1 −σ1 0 ] , . . . , [ 0 σK/2 −σK/2 0 ]) ,
X̂ := diag ( ρ1, . . . , ρK , [ σ1 0
0 σ1
] , . . . , [ σK/2 0
0 σK/2
]) ,
so that L = ZXZ> and L̂ = ZX̂Z>. Applying the Cauchy-Binet formula twice, we can write the determinant of the principal submatrices of both L and L̂:
det(LY ) = ∑
I∈([2K]k ) ∑ J∈([2K]k ) det(XI,J) det(ZY,I) det(ZY,J), (18)
det(L̂Y ) = ∑
I∈([2K]k ) ∑ J∈([2K]k ) det(X̂I,J) det(ZY,I) det(ZY,J)
= ∑
I∈([2K]k )
det(X̂I) det(ZY,I) 2, (19)
where Eq. (19) follows from the fact that X̂ is diagonal, which means that det(X̂I,J) = 0 for I 6= J . When the size of Y is equal to the rank of L (i.e., k = 2K), the summations in Eqs. (18) and (19) simplify to single terms: det(LY ) = det(X) det(ZY,:)2 and det(L̂Y ) = det(X̂) det(ZY,:)2. Now, observe that the determinants of the full X and X̂ matrices are identical: det(X) = det(X̂) =∏K i=1 ρi ∏K/2 j=1 σ 2 j . Hence, it holds that det(LY ) = det(L̂Y ). This proves the second statement of the theorem.
To prove that det(LY ) ≤ det(L̂Y ) for smaller subsets Y , we will use the following: Claim 1. For every I, J ∈ ( [2K] k ) such that det(XI,J) 6= 0, there exists a (nonempty) collection of
subset pairs S(I, J) ⊆ ( [2K] k ) × ( [2K] k ) such that∑
(I′,J′)∈S(I,J)
det(XI,J) det(ZY,I) det(ZY,J) ≤ ∑
(I′,J′)∈S(I,J)
det(X̂I,I) det(ZY,I) 2. (20)
Claim 2. The number of nonzero terms in Eq. (18) is identical to that in Eq. (19).
Combining Claim 1 with Claim 2 yields det(LY ) = ∑
I,J∈([2K]k )
det(XI,J) det(ZY,I) det(ZY,J) ≤ ∑
I∈([2K]k )
det(X̂I,I) det(ZY,I) 2 = det(L̂Y ).
We conclude the proof of Theorem 1. Below we provide proofs for Claim 1 and Claim 2.
Proof of Claim 1. Recall that X is a block-diagonal matrix, where each block is of size either 1-by-1, containing ρi, or 2-by-2, containing both σj and −σj in the form [ 0 σj −σj 0 ] . A submatrix XI,J ∈ Rk×k with rows I and columns J will only have a nonzero determinant if it contains no
all-zero row or column. Hence, any XI,J with nonzero determinant will have the following form (or some permutation of this block-diagonal):
XI,J =  ρp1 · · · 0 ... . . . ... 0 0 . . . ρp|PI,J | ±σq1 · · · 0 ... . . . ... 0 · · · ±σq|QI,J | 0 σr1 −σr1 0 . . .
0 0 σr|RI,J | −σr|RI,J | 0  (21)
and we denote P I,J := {p1, . . . , p|P I,J |}, QI,J := {q1, . . . , q|QI,J |}, and RI,J := {r1, . . . , r|RI,J |}. Indices p ∈ P I,J yield a diagonal matrix with entries ρp. For such p, both I and J must contain index p. Indices r ∈ RI,J yield a block-diagonal matrix of the form [ 0 σr −σr 0 ] . For such r, both I and J must contain a pair of indices, (K + 2r − 1,K + 2r). Finally, indices q ∈ QI,J yield a diagonal matrix with entries of ±σq (the sign can be + or −). For such q, I contains K + 2q − 1 or K + 2q, and J must contain the other. Note that there is no intersection between QI,J and RI,J .
If QI,J is an empty set (i.e., I = J), then det(XI,J) = det(X̂I,J) and
det(XI,J) det(ZY,I) det(ZY,J) = det(X̂I) det(ZY,I) 2. (22)
Thus, the terms in Eq. (18) in this case appear in Eq. (19). Now assume that QI,J 6= ∅ and consider the following set of pairs:
S(I, J) := {(I ′, J ′) : P I,J = P I′,J′ , QI,J = QI′,J′ , RI,J = RI′,J′}.
In other words, for (I ′, J ′) ∈ S(I, J), the diagonal XI′,J′ contains ρp, [ 0 σr −σr 0 ] exactly as in XI,J . However, the signs of the σr’s may differ from XI,J . Combining this observation with the definition of X̂ ,
|det(XI′,J′)| = |det(XI,J)| = det(X̂I) = det(X̂I′) = det(X̂J) = det(X̂J′). (23) Therefore, ∑
(I′,J′)∈S(I,J)
det(XI′,J′) det(ZY,I′) det(ZY,J′) (24)
≤ ∑
(I′,J′)∈S(I,J)
|det(XI′,J′)|det(ZY,I′) det(ZY,J′) (25)
= det(X̂I) ∑
(I′,J′)∈S(I,J)
det(ZY,I′) det(ZY,J′) (26)
≤ det(X̂I) ∑
(I′,∗)∈S(I,J)
det(ZY,I′) 2 (27)
= ∑
(I′,∗)∈S(I,J)
det(X̂I′) det(ZY,I′) 2 (28)
where the third line comes from Eq. (23) and the fourth line follows from the rearrangement inequality. Note that application of this inequality does not change the number of terms in the sum. This completes the proof of Claim 1.
Proof of Claim 2. In Eq. (19), observe that det(X̂I) det(ZY,I)2 6= 0 if and only if det(X̂I) 6= 0. Since all ρi’s and σj’s are positive, the number of I ⊆ [2K], |I| = k such that det(X̂I) 6= 0 is equal to ( 2K k ) . Similarly, the number of nonzero terms in Eq. (18) equals the number of possible
choices of I, J ∈ ( [2K] k ) such that det(XI,J) 6= 0. This can be counted as follows: first choose i
items in {ρ1, . . . , ρK} for i = 0, . . . , k; then, choose j items in {[ 0 σ1 −σ1 0 ] , . . . , [ 0 σK/2 −σK/2 0 ]} for j = 0, . . . , bk−i2 c; lastly, choose k − i − 2j of {±σq; q /∈ RI,J}, then choose the sign for each of these (σq or −σq). Combining all of these choices, the total number of nonzero terms is:
k∑ i=0 ( K i ) ︸ ︷︷ ︸
choice of ρp
b k−i2 c∑ j=0
( K/2
j ) ︸ ︷︷ ︸
choice of [
0 σr −σr 0
] ( K/2− j k − i− 2j ) 2k−i−2j︸ ︷︷ ︸
choice of±σq
(29)
= k∑ i=0 ( K i ) ( K k − i ) (30)
=
( 2K
k
) (31)
where the second line comes from the fact that ( 2n m ) = ∑bm2 c j=0 ( n j )( n−j m−2j ) 2m−2j for any integers n,m ∈ N such thatm ≤ 2n (see (1.69) in Quaintance (2010)), and the third line follows from the fact that ∑r i=0 ( m i )( n r−i ) = ( n+m r ) for n,m, r ∈ N (Vandermonde’s identity). Hence, both the number
of nonzero terms in Eqs. (18) and (19) is equal to ( 2K k ) . This completes the proof of Claim 2.
E.2 PROOF OF PROPOSITION 1
Proposition 1. The tree-based sampling procedure SAMPLEDPP in Algorithm 3 runs in time O(K + k3 logM + k4), where k is the size of the sampled set***.
Proof of Proposition 1. Since computing p` takes O(k2) from Eq. (12), and since the binary tree has depth O(logM), SAMPLEITEM in Algorithm 3 runs in O(k2 logM) time. Moreover, the query matrix QY can be updated in O(k3) time as it only requires a k-by-k matrix inversion. Therefore, the overall runtime of the tree-based elementary DPP sampling algorithm (after pre-processing) is O(k3 logM + k4). This improves the previous O(k4 logM) runtime studied in Gillenwater et al. (2019). Combining this with elementary DPP selection (Line 15 in Algorithm 3), we can sample a set in O(K + k3 logM + k4) time. This completes the proof of Proposition 1.
E.3 PROOF OF THEOREM 2
Theorem 2. Given an NDPP kernel L = V V > + B(D −D>)B> for V ,B ∈ RM×K ,D ∈ RK×K , consider the proposal kernel L̂ as proposed in Section 4.1. Let {σj}K/2j=1 be the positive eigenvalues obtained from the Youla decomposition of B(D−D>)B>. If V ⊥ B, then det(L̂+I)det(L+I) =∏K/2 j=1 ( 1 +
2σj σ2j+1 ) ≤ (1 + ω)K/2, where ω = 2K ∑K/2 j=1 2σj σ2j+1 ∈ (0, 1].
Proof of Theorem 2. Since the column spaces of V and B are orthogonal, the corresponding eigenvectors are also orthogonal, i.e., Z>Z = I2K . Then,
det(L+ I) = det(ZXZ> + I) = det(XZ>Z + I2K) = det(X + I2K) (32)
= K∏ i=1 (ρi + 1) K/2∏ j=1 det ([ 1 σj −σj 1 ]) (33)
= K∏ i=1 (ρi + 1) K/2∏ j=1 (σ2j + 1) (34)
***Computing p` via Eq. (12) improves on Gillenwater et al. (2019)’s O(k4 logM) runtime for this step.
and similarly
det(L̂+ I) = K∏ i=1 (ρi + 1) K/2∏ j=1 (σj + 1) 2. (35)
Combining Eqs. (34) and (35), we have that
det(L̂+ I) det(L+ I) = K/2∏ j=1 (σj + 1) 2 (σ2j + 1) = K/2∏ j=1
( 1 +
2σj σ2j + 1
) ≤ 1 + 2 K K/2∑ j=1 2σj σ2j + 1 K/2 (36) where the inequality holds from the Jensen’s inequality. This completes the proof of Theorem 2.

Review summaries:
This paper studies scalable sampling methods for NDPPs: - By assuming low-rank structure of the kernel, the paper proposes a linear-time sampling method, which is faster than previous sampling algorithm with cubic runtime for general kernel. - Furthermore, this paper develops a scalable sub-linear-time rejection sampling algorithm for a subclass of low-rank NDPPs that are called ONDPPs. - Through experiments, it is shown that the predictive performance of ONDPPs is not degraded compared to NDPPs. By adding a regularization term in the optimization, the rejection probability is greatly reduced.

This paper studies the problem of sampling from non-symmetrical Determinantal Point Processes (NDPPs); in particular, this paper focuses on exact sampling for low-rank NDPPs.  The main contributions of this paper are: 1. This paper proposes to adapt the Cholesky-decomposition-based sampler for DPPs to a linear-time sampler for low-rank NDPPs (rank of NDPP << size of the ground set) 2. This paper proposes to use rejection sampling to implement a "sub-linear-time" (assuming that the # of rejections are bounded by a small constant) sampler for a subclass of NDPPs called orthogonal NDPPs (ONDPPs), by leveraging an existing sub-linear-time sampler for DPPs. 3. This paper shows that empirically, in terms of modeling real-world datasets, ONDPPs are as effective as the general NDPPs and more importantly, we can efficiently learn ONDPPs in a way such that when it comes to rejection sampling, the # of rejections is bounded by a small number.

This paper gives an exact sampling algorithm to sample from DPPs based on non-symmetrical, low-rank, kernel matrices of the form $L=VV^\top + B(D-D^\top)B^\top$ where $V$ and $B$ are of size $M$ by $K$ and $D$ is square, non-symmetric and of size $K$ by $K$. The marginal kernel of such a DPP may be written in the form $K=ZWZ^\top$ with $Z$ a $M$ by $2K$ matrix and $W$ a square matrix of size $2K$ that is not symmetric. The authors then: - discuss in Section 2 how Poulson's algorithm is adapted to sample from such NDPPS - in Section 3: show that, as L is low-rank, the generic $O(M^3)$ cost of Poulson's algorithm can be reduced to $O(MK^2$ as updates of the marginal kernel can be efficiently done by only updating the inner $W$ matrix - in Section 4: give a sublinear-time rejection-sampling based algorithm adapting the tree-based algorithm of Gillenwater et al. In this section, they: 	* propose a well-adapted, DPP based on a symmetric kernel $\hat{L}$ that is easy-to-sample (it is a DPP based on a symmetric PSD kernel and can be sampled from many different fast existing algorithms), and for which any set $Y$ verifies $det(L_Y)\leq det(\hat{L}_Y)$ where the upper-bound is actually attained (Thm 1). In order to exactly control the rejection rate of their proposed algorithm (that is equal to $det(\hat{L}+I) / det(L+I)$, the authors simplify the model by supposing that $V$ and $B$ are orthogonal to each other (Thm 2) 	* in passing, the authors observe that the complexity of Gillenwater et al's algorithm can in fact be easily improved by a factor $k$ where $k$ is the size of the sampled set  - in Sections 5 and 6, experiments are provided comparing their samping method with the state-of-the-art