# Kronecker Determinantal Point Processes

Zelda Mariet
Massachusetts Institute of Technology
Cambridge, MA 02139
<EMAIL_ADDRESS>
Suvrit Sra
Massachusetts Institute of Technology
Cambridge, MA 02139
<EMAIL_ADDRESS>

###### Abstract

Determinantal Point Processes (DPPs) are probabilistic models over all subsets
a ground set of $N$ items. They have recently gained prominence in several
applications that rely on “diverse” subsets. However, their applicability to
large problems is still limited due to the $\mathcal{O}(N^{3})$ complexity of
core tasks such as sampling and learning. We enable efficient sampling and
learning for DPPs by introducing KronDpp, a DPP model whose kernel matrix
decomposes as a tensor product of multiple smaller kernel matrices. This
decomposition immediately enables fast _exact_ sampling. But contrary to what
one may expect, leveraging the Kronecker product structure for speeding up DPP
learning turns out to be more difficult. We overcome this challenge, and
derive batch and stochastic optimization algorithms for efficiently learning
the parameters of a KronDpp.

##  1 Introduction

Determinantal Point Processes (DPPs) are discrete probability models over the
subsets of a ground set of $N$ items. They provide an elegant model to assign
probabilities to an exponentially large sample, while permitting tractable
(polynomial time) sampling and marginalization. They are often used to provide
models that balance “diversity” and quality, characteristics valuable to
numerous problems in machine learning and related areas [17].

The antecedents of DPPs lie in statistical mechanics [24], but since the
seminal work of [15] they have made inroads into machine learning. By now they
have been applied to a variety of problems such as document and video
summarization [21, 6], sensor placement [14], recommender systems [31], and
object retrieval [2]. More recently, they have been used to compress fully-
connected layers in neural networks [26] and to provide optimal sampling
procedures for the Nyström method [20]. The more general study of DPP
properties has also garnered a significant amount of interest, see e.g., [12,
16, 17, 1, 7, 18, 5, 23].

However, despite their elegance and tractability, widespread adoption of DPPs
is impeded by the $O(N^{3})$ cost of basic tasks such as (exact) sampling [12,
17] and learning [12, 17, 10, 25]. This cost has motivated a string of recent
works on approximate sampling methods such as MCMC samplers [13, 20] or core-
set based samplers [19]. The task of learning a DPP from data has received
less attention; the methods of [10, 25] cost $O(N^{3})$ per iteration, which
is clearly unacceptable for realistic settings. This burden is partially
ameliorated in [9], who restrict to learning low-rank DPPs, though at the
expense of being unable to sample subsets larger than the chosen rank.

These considerations motivate us to introduce KronDpp, a DPP model that uses
Kronecker (tensor) product kernels. As a result, KronDpp enables us to learn
large sized DPP kernels, while also permitting efficient (exact and
approximate) sampling. The use of Kronecker products to scale matrix models is
a popular and effective idea in several machine-learning settings [28, 27, 8,
30]. But as we will see, its efficient execution for DPPs turns out to be
surprisingly challenging.

To make our discussion more concrete, we recall some basic facts now. Suppose
we have a ground set of $N$ items $\mathcal{Y}=\\{1,\ldots,N\\}$. A discrete
DPP over $\mathcal{Y}$ is a probability measure $\mathcal{P}$ on
$2^{\mathcal{Y}}$ parametrized by a positive definite matrix $K$ (the
_marginal kernel_) such that $0\preceq K\preceq I$, so that for any
$Y\in\mathcal{Y}$ drawn from $\mathcal{P}$, the measure satisfies

$\forall A\subseteq\mathcal{Y},\qquad\mathcal{P}(A\subseteq Y)=\det(K_{A}),$
(1)

where $K_{A}$ is the submatrix of $K$ indexed by elements in $A$ (i.e.,
$K_{A}=[K_{ij}]_{i,j\in A}$). If a DPP with marginal kernel $K$ assigns
nonzero probability to the empty set, the DPP can alternatively be
parametrized by a positive definite matrix $L$ (the _DPP kernel_) so that

$\mathcal{P}(Y)\propto\det(L_{Y})\quad\implies\quad\mathcal{P}(Y)=\frac{\det(L_{Y})}{\det(L+I)}.$
(2)

A brief manipulation (see e.g., [17, Eq. 15]) shows that when the inverse
exists, $L=K(I-K)^{-1}$. The determinants, such as in the normalization
constant in (2), make operations over DPPs typically cost
$\mathcal{O}(N^{3})$, which is a key impediment to their scalability.

Therefore, if we consider a class of DPP kernels whose structure makes it easy
to compute determinants, we should be able to scale up DPPs. An alternative
approach towards scalability is to restrict the size of the subsets, as done
in $k$-DPP [16] or when using rank-$k$ DPP kernels [9] (where $k\ll N$). Both
these approaches still require $O(N^{3})$ preprocessing for exact sampling;
another caveat is that they limit the DPP model by assigning zero
probabilities to sets of cardinality greater than $k$.

In contrast, KronDpp uses a kernel matrix of the form
$L=L_{1}\otimes\ldots\otimes L_{m}$, where each _sub-kernel_ $L_{i}$ is a
smaller positive definite matrix. This decomposition has two key advantages:
(i) it significantly lowers the number of parameters required to specify the
DPP from $N^{2}$ to $\mathcal{O}(N^{2/m})$ (assuming the sub-kernels are
roughly the same size); and (ii) it enables fast sampling and learning.

For ease of exposition, we describe specific details of KronDpp for $m=2$; as
will become clear from the analysis, typically the special cases $m=2$ and
$m=3$ should suffice to obtain low-complexity sampling and learning
algorithms.

##### Contributions.

Our main contribution is the KronDpp model along with efficient algorithms for
sampling from it and learning a Kronecker factored kernel. Specifically,
inspired by the algorithm of [25], we develop KrK-Picard (Kronecker-Kernel
Picard), a block-coordinate ascent procedure that generates a sequence of
Kronecker factored estimates of the DPP kernel while ensuring monotonic
progress on its (difficult, nonconvex) objective function. More importantly,
we show how to implement KrK-Picard to run in $\mathcal{O}(N^{2})$ time when
implemented as a batch method, and in $\mathcal{O}(N^{3/2})$ time and
$\mathcal{O}(N)$ space, when implemented as a stochastic method. As alluded to
above, unlike many other uses of Kronecker models, KronDpp does not admit
trivial scaling up, largely due to extensive dependence of DPPs on arbitrary
submatrices of the DPP kernel. An interesting theoretical nugget that arises
from our analysis is the combinatorial problem that we call _subset
clustering_ , a problem whose (even approximate) solution can lead to further
speedups of our algorithms.

##  2 Preliminaries

We begin by recalling basic properties of Kronecker products needed in our
analysis; we omit proofs of these well-known results for brevity. The
Kronecker (tensor) product of $A\in\mathbb{R}^{p\times q}$ with
$B\in\mathbb{R}^{r\times s}$ two matrices is defined as the $pr\times qs$
block matrix $A\otimes B=[a_{ij}B]_{i,j=1}^{p,q}$.

We denote the block $a_{ij}B$ in $A\otimes B$ by $(A\otimes B)_{(ij)}$ for any
valid pair $(i,j)$, and extend the notation to non-Kronecker product matrices
to indicate the submatrix of size $r\times s$ at position $(i,j)$.

######  Proposition 2.1.

Let $A,B,C,D$ be matrices of sizes so that $AC$ and $BD$ are well-defined.
Then,

  1. (i)

If $A,B\succeq 0$, then, $A\otimes B\succeq 0$;

  2. (ii)

If $A$ and $B$ are invertible then so is $A\otimes B$, with $(A\otimes
B)^{-1}=A^{-1}\otimes B^{-1}$;

  3. (iii)

$(A\otimes B)(C\otimes D)$ = $(AC)\otimes(BD)$.

An important consequence of Prop. 2.1$(iii)$ is the following corollary.

######  Corollary 2.2.

Let $A=P_{A}D_{A}P_{A}^{\top}$ and $B=P_{B}D_{B}P_{B}^{\top}$ be the
eigenvector decompositions of $A$ and $B$. Then, $A\otimes B$ diagonalizes as
$(P_{A}\otimes P_{B})(D_{A}\otimes D_{B})(P_{A}\otimes P_{B})^{\top}$.

We will also need the notion of partial trace operators, which are perhaps
less well-known:

######  Definition 2.3.

Let $A\in\mathbb{R}^{N_{1}N_{2}\times N_{1}N_{2}}$. The _partial traces_
$\operatorname{Tr}_{1}(A)$ and $\operatorname{Tr}_{2}(A)$ are defined as
follows:

$\operatorname{Tr}_{1}(A):=\left[\operatorname{Tr}(A_{(ij)}\right]_{1\leq
i,j\leq N_{1}}\in\mathbb{R}^{N_{1}\times
N_{1}},\qquad\operatorname{Tr}_{2}(A):=\sum\nolimits_{i=1}^{N_{1}}A_{(ii)}\in\mathbb{R}^{N_{2}\times
N_{2}}.$

The action of partial traces is easy to visualize: indeed,
$\operatorname{Tr}_{1}(A\otimes B)=\operatorname{Tr}(B)A$ and
$\operatorname{Tr}_{2}(A\otimes B)=\operatorname{Tr}(A)B$. For us, the most
important property of partial trace operators is their positivity.

######  Proposition 2.4.

$\operatorname{Tr}_{1}$ and $\operatorname{Tr}_{2}$ are positive operators,
i.e., for $A\succ 0$, $\operatorname{Tr}_{1}(A)\succ 0$ and
$\operatorname{Tr}_{2}(A)\succ 0$.

###### Proof.

Please refer to [4, Chap. 4]. ∎

##  3 Learning the kernel matrix for KronDpp

In this section, we consider the key difficult task for KronDpps: learning a
Kronecker product kernel matrix from $n$ observed subsets
$Y_{1},\ldots,Y_{n}$. Using the definition (2) of $\mathcal{P}(Y_{i})$,
maximum-likelihood learning of a DPP with kernel $L$ results in the
optimization problem:

$\arg\max_{L\succ
0}\quad\phi(L),\qquad\phi(L)=\frac{1}{n}\sum_{i=1}^{n}\left(\log\det(L_{Y_{i}})-\log\det(L+I)\right).$
(3)

This problem is nonconvex and conjectured to be NP-hard [15, Conjecture 4.1].
Moreover the constraint $L\succ 0$ is nontrivial to handle. Writing $U_{i}$ as
the indicator matrix for $Y_{i}$ of size $N\times|Y_{i}|$ so that
$L_{Y_{i}}=U_{i}^{\top}LU_{i}$, the gradient of $\phi$ is easily seen to be

$\Delta:=\nabla\phi(L)=\frac{1}{n}\sum\nolimits_{i=1}^{n}U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}-(L+I)^{-1}.$
(4)

In [25], the authors derived an iterative method (“the Picard iteration”) for
computing an $L$ that solves $\Delta=0$ by running the simple iteration

$L\leftarrow L+L\Delta L.$ (5)

Moreover, iteration (5) is guaranteed to monotonically increase the log-
likelihood $\phi$ [25]. But these benefits accrue at a cost of $O(N^{3})$ per
iteration, and furthermore a direct application of (5) cannot guarantee the
Kronecker structure required by KronDpp.

###  3.1 Optimization algorithm

Our aim is to obtain an efficient algorithm to (locally) optimize (3). Beyond
its nonconvexity, the Kronecker structure $L=L_{1}\otimes L_{2}$ imposes
another constraint. As in [25] we first rewrite $\phi$ as a function of
$S=L^{-1}$, and re-arrange terms to write it as

$\phi(S)=\underbrace{\vphantom{\sum\nolimits_{a}^{b}\frac{1}{n}}\log\det(S)}_{f(S)}+\underbrace{\frac{1}{n}\sum\nolimits_{i=1}^{n}\log\det\left(U_{i}^{\top}S^{-1}U_{i}\right)-\log\det(I+S)}_{g(S)}.$
(6)

It is easy to see that $f$ is concave, while a short argument shows that $g$
is convex [25]. An appeal to the convex-concave procedure [29] then shows that
updating $S$ by solving $\nabla\\!f(S^{(k+1)})+\nabla\\!g(S^{(k)})=0$, which
is what (5) does [25, Thm. 2.2], is guaranteed to monotonically increase
$\phi$.

But for KronDpp this idea does not apply so easily: due the constraint
$L=L_{1}\otimes L_{2}$ the function

$g_{\tiny\otimes}:(S_{1},S_{2})\rightarrow\tfrac{1}{n}\sum\nolimits_{i=1}^{n}\log\det\left(U_{i}^{\top}(S_{1}\otimes
S_{2})^{-1}U_{i}\right)-\log\det(I+S_{1}\otimes S_{2}),$

fails to be convex, precluding an easy generalization. Nevertheless, for fixed
$S_{1}$ or $S_{2}$ the functions

$\begin{cases}f_{1}:S_{1}\mapsto f(S_{1}\otimes S_{2})\\\ g_{1}:S_{1}\mapsto
g(S_{1}\otimes S_{2})\\\ \end{cases},\qquad\begin{cases}f_{2}:S_{2}\rightarrow
f(S_{1}\otimes S_{2})\\\ g_{2}:S_{2}\rightarrow g(S_{1}\otimes S_{2})\\\
\end{cases}$

are once again concave or convex. Indeed, the map $\otimes:S_{1}\rightarrow
S_{1}\otimes S_{2}$ is linear and $f$ is concave, and $f_{1}=f\circ\otimes$ is
also concave; similarly, $f_{2}$ is seen to be concave and $g_{1}$ and $g_{2}$
are convex. Hence, by generalizing the arguments of [29, Thm. 2] to our
“block-coordinate” setting, updating via

$\nabla\\!f_{i}\bigl{(}{S_{i}}^{(k+1)}\bigr{)}=-\nabla\\!g_{i}\bigl{(}{S_{i}}^{(k)}\bigr{)},\quad\text{for}\
i=1,2,$ (7)

should increase the log-likelihood $\phi$ at each iteration. We prove below
that this is indeed the case, and that updating as per (7) ensure positive
definiteness of the iterates as well as monotonic ascent.

####  3.1.1 Positive definite iterates and ascent

In order to show the positive definiteness of the solutions to (7), we first
derive their closed form.

######  Proposition 3.1 (Positive definite iterates).

For $S_{1}\succ 0$, $S_{2}\succ 0$, the solutions to (7) are given by the
following expressions:

$\displaystyle\nabla\\!f_{1}(X)=-\nabla\\!g_{1}(S_{1})$ $\displaystyle\iff
X^{-1}=\operatorname{Tr}_{1}\\!\left((I\otimes S_{2})(L+L\Delta
L)\right)/N_{2}$ $\displaystyle\nabla\\!f_{2}(X)=-\nabla\\!g_{2}(S_{2})$
$\displaystyle\iff X^{-1}=\operatorname{Tr}_{2}\left((S_{1}\otimes
I)(L+L\Delta L)\right)/N_{1}.$

Moreover, these solutions are positive definite.

###### Proof.

The details are somewhat technical, and are hence given in Appendix A. We know
that $L\succ 0\implies L+L\Delta L\geq 0$, because $L-L(I+L)^{-1}L\succ 0$.
Since the partial trace operators are positive (Prop. 2.4), it follows that
the solutions to (7) are also positive definite. ∎

We are now ready to establish that these updates ensure monotonic ascent in
the log-likelihood.

######  Theorem 3.2 (Ascent).

Starting with $L_{1}^{(0)}\succ 0$, $L_{2}^{(0)}\succ 0$, updating according
to (7) generates positive definite iterates $L_{1}^{(k)}$ and $L_{2}^{(k)}$,
and the sequence $\bigl{\\{}\phi\bigl{(}L_{1}^{(k)}\otimes
L_{2}^{(k)}\bigr{)}\bigr{\\}}_{k\geq 0}$ is non-decreasing.

###### Proof.

Updating according to (7) generates positive definite matrices $S_{i}$, and
hence positive definite subkernels $L_{i}=S_{i}$. Moreover, due to the
convexity of $g_{1}$ and concavity of $f_{1}$, for matrices $A,B\succ 0$

$\displaystyle f_{1}(B)\leq f_{1}(A)+\nabla\\!f_{1}(A)^{\top}(B-A),$
$\displaystyle g_{1}(A)\geq g_{1}(B)+\nabla\\!g_{1}(B)^{\top}(A-B).$

Hence, $f_{1}(A)+g_{1}(A)\geq
f_{1}(B)+g_{1}(B)+(\nabla\\!f_{1}(A)+\nabla\\!g_{1}(B))^{\top}(A-B)$.

Thus, if $S_{1}^{(k)},S_{1}^{(k+1)}$ verify (7), by setting $A=S_{1}^{(k+1)}$
and $B=S_{1}^{(k)}$ we have

$\phi\bigl{(}L_{1}^{(k+1)}\otimes
L_{2}^{(k)}\bigr{)}=f_{1}\\!\bigl{(}S_{1}^{(k+1)}\bigr{)}+g_{1}\\!\bigl{(}S_{1}^{(k+1)}\bigr{)}\geq
f_{1}\\!\bigl{(}S_{1}^{(k)}\bigr{)}+g_{1}\\!\bigl{(}S_{1}^{(k)}\bigr{)}=\phi\bigl{(}L_{1}^{(k)}\otimes
L_{2}^{(k)}\bigr{)}.$

The same reasoning holds for $L_{2}$, which proves the theorem. ∎

As $\operatorname{Tr}_{1}\\!\left((I\otimes S_{2})L\right)=N_{2}L_{1}$ (and
similarly for $L_{2}$), updating as in (7) is equivalent to updating

$L_{1}\leftarrow L_{1}+\operatorname{Tr}_{1}\\!\left((I\otimes
L^{-1}_{2})(L\Delta L)\right)/N_{2},\qquad L_{2}\leftarrow
L_{2}+\operatorname{Tr}_{2}\left((L^{-1}_{1}\otimes I)(L\Delta
L)\right)/N_{1}.$

Genearlization. We can generalize the updates to take an additional step-size
parameter $a$:

$L_{1}\leftarrow L_{1}+a\operatorname{Tr}_{1}\\!\left((I\otimes
L^{-1}_{2})(L\Delta L)\right)/N_{2},\qquad L_{2}\leftarrow
L_{2}+a\operatorname{Tr}_{2}\left((L^{-1}_{1}\otimes I)(L\Delta
L)\right)/N_{1}.$

Experimentally, $a>1$ (as long as the updates remain positive definite) can
provide faster convergence, although the monotonicity of the log-likelihood is
no longer guaranteed. We found experimentally that the range of admissible $a$
is larger than for Picard, but decreases as $N$ grows larger.

The arguments above easily generalize to the multiblock case. Thus, when
learning $L=L_{1}\otimes\cdots\otimes L_{m}$, by writing $E_{ij}$ the matrix
with a 1 in position $(i,j)$ and zeros elsewhere, we update $L_{k}$ as

$(L_{k})_{ij}\leftarrow(L_{k})_{ij}+N_{k}/(N_{1}\ldots
N_{m})\operatorname{Tr}\left[(L_{1}\otimes\ldots\otimes L_{k-1}\otimes
E_{ij}\otimes L_{k+1}\otimes\ldots\otimes L_{m})(L\Delta L)\right].$

From the above updates it is not transparent whether the Kronecker product
saves us any computation. In particular, it is not clear whether the updates
can be implemented to run faster than $O(N^{3})$. We show below in the next
section how to implement these updates efficiently.

####  3.1.2 Algorithm and complexity analysis

From Theorem 3.2, we obtain Algorithm 1 (which is different from the Picard
iteration in [25], because it operates alternatingly on each subkernel). It is
important to note that a further speedup to Algorithm 1 can be obtained by
performing stochastic updates, i.e., instead of computing the full gradient of
the log-likelihood, we perform our updates using only one (or a small
minibatch) subset $Y_{i}$ at each step instead of iterating over the entire
training set; this uses the stochastic gradient
$\Delta=U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}-(I+L)^{-1}$.

Algorithm 1 KrK-Picard iteration

Input: Matrices $L_{1},L_{2}$, training set $T$, parameter $a$.

for $i=1$ to maxIter do

$L_{1}\leftarrow L_{1}+a\operatorname{Tr}_{1}\\!\left((I\otimes
L^{-1}_{2})(L\Delta L)\right)/N_{2}$ // or update stochastically

$L_{2}\leftarrow L_{2}+a\operatorname{Tr}_{2}\left((L^{-1}_{1}\otimes
I)(L\Delta L)\right)/N_{1}$ // or update stochastically

end forreturn $(L_{1},L_{2})$

The crucial strength of Algorithm 1 lies in the following result:

######  Theorem 3.3 (Complexity).

For $N_{1}\approx N_{2}\approx\sqrt{N}$, the updates in Algorithm 1 can be
computed in $\mathcal{O}(n\kappa^{3}+N^{2})$ time and $\mathcal{O}(N^{2})$
space, where $\kappa$ is the size of the largest training subset. Furthermore,
stochastic updates can be computed in $\mathcal{O}(N\kappa^{2}+N^{3/2})$ time
and $\mathcal{O}(N+\kappa^{2})$ space.

Indeed, by leveraging the properties of the Kronecker product, the updates can
be obtained without computing $L\Delta L$. This result is non-trivial: the
components of $\Delta$, $\frac{1}{n}\sum_{i}U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}$
and $(I+L)^{-1}$, must be considered separately for computational efficiency.
The proof is provided in App. B. However, it seems that considering more than
2 subkernels does not lead to further speed-ups.

If $N_{1}\approx N_{2}\approx\sqrt{N}$, these complexities become:

  * –

for non-stochastic updates: $\mathcal{O}(n\kappa^{3}+N^{2})$ time,
$\mathcal{O}(N^{2})$ space,

  * –

for stochastic updates: $\mathcal{O}(N\kappa^{3}+N^{3/2})$ time,
$\mathcal{O}(\kappa^{2}+N)$ space.

This is a marked improvement over [25], which runs in $\mathcal{O}(N^{2})$
space and $\mathcal{O}(n\kappa^{3}+N^{3})$ time (non-stochastic) or
$\mathcal{O}(N^{3})$ time (stochastic); Algorithm 1 also provides faster
stochastic updates than [9]. However, one may wonder if by learning the sub-
kernels by alternating updates the log-likelihood converges to a sub-optimal
limit. The next section discusses how to jointly update $L_{1}$ and $L_{2}$.

###  3.2 Joint updates

We also analyzed the possibility of updating $L_{1}$ and $L_{2}$ jointly: we
update $L\leftarrow L+L\Delta L$ and then recover the Kronecker structure of
the kernel by defining the updates $L_{1}^{\prime}$ and $L_{2}^{\prime}$ such
that:

$\begin{cases}(L_{1}^{\prime},L_{2}^{\prime})\text{ minimizes }\|L+L\Delta
L-L_{1}^{\prime}\otimes L_{2}^{\prime}\|^{2}_{F}\\\ L_{1}^{\prime}\succ
0,L_{2}^{\prime}\succ 0,\|L_{1}^{\prime}\|=\|L_{2}^{\prime}\|\end{cases}$ (8)

We show in appendix C that such solutions exist and can be computed by from
the first singular value and vectors of the matrix
$R=\left[\operatorname{vec}((L^{-1}+\Delta)_{(ij)})^{\top}\right]_{i,j=1}^{N_{1}}$.
Note however that in this case, there is no guaranteed increase in log-
likelihood. The pseudocode for the related algorithm (Joint-Picard) is given
in appendix C.1. An analysis similar to the proof of Thm. 3.3 shows that the
updates can be obtained $\mathcal{O}(n\kappa^{3}+\max(N_{1},N_{2})^{4})$.

###  3.3 Memory-time trade-off

Although KronDpps have tractable learning algorithms, the memory requirements
remain high for non-stochastic updates, as the matrix
$\Theta=\frac{1}{n}\sum_{i}U_{i}L^{-1}_{Y_{i}}U_{i}^{\top}$ needs to be
stored, requiring $\mathcal{O}(N^{2})$ memory. However, if the training set
can be subdivised such that

$\\{Y_{1},\ldots,Y_{n}\\}=\cup_{k=1}^{m}S_{k}\quad s.t.\quad\forall
k,\left|\cup_{Y\in S_{k}}Y\right|<z,$ (9)

$\Theta$ can be decomposed as $\frac{1}{n}\sum_{k=1}^{m}\Theta_{k}$ with
$\Theta_{k}=\sum_{Y_{i}\in S_{k}}U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}$. Due to the
bound in Eq. 9, each $\Theta_{k}$ will be sparse, with only $z^{2}$ non-zero
coefficients. We can then store each $\Theta_{k}$ with minimal storage and
update $L_{1}$ and $L_{2}$ in $\mathcal{O}(n\kappa^{3}+mz^{2}+N^{3/2})$ time
and $\mathcal{O}(mz^{2}+N)$ space.

Determining the existence of such a partition of size $m$ is a variant of the
NP-Hard Subset-Union Knapsack Problem (SUKP) [11] with $m$ knapsacks and where
the value of each item (i.e. each $Y_{i}$) is equal to 1: a solution to SUKP
of value $n$ with $m$ knapsacks is equivalent to a solution to 9. However, an
approximate partition can also be simply constructed via a greedy algorithm.

##  4 Sampling

Sampling exactly (see Alg. 2 and [17]) from a full DPP kernel costs
$\mathcal{O}(N^{3}+Nk^{3})$ where $k$ is the size of the sampled subset. The
bulk of the computation lies in the initial eigendecomposition of $L$; the $k$
orthonormalizations cost $\mathcal{O}(Nk^{3})$. Although the
eigendecomposition need only happen once for many iterations of sampling,
exact sampling is nonetheless intractable in practice for large $N$.

Algorithm 2 Sampling from a DPP kernel $L$

Input: Matrix $L$.

Eigendecompose $L$ as $\\{(\lambda_{i},v_{i})\\}_{1\leq i\leq N}$.

$J\leftarrow\emptyset$

for $i=1$ to $N$ do

$J\rightarrow J\cup\\{i\\}$ with probability $\lambda_{i}/(\lambda_{i}+1)$.

end for

$V\leftarrow\\{v_{i}\\}_{i\in J}$, $Y\leftarrow\emptyset$

while $|V|>0$ do

Sample $i$ from $\\{1\ldots N\\}$ with probability $\frac{1}{|V|}\sum_{v\in
V}v_{i}^{2}$

$Y\leftarrow Y\cup\\{i\\}$, $V\leftarrow V_{\bot}$, where $V_{\bot}$ is an
orthonormal basis of the subspace of $V$ orthonormal to $e_{i}$

end whilereturn $Y$

It follows from Prop. 2.2 that for KronDpps, the eigenvalues $\lambda_{i}$ can
be obtained in $\mathcal{O}(N_{1}^{3}+N_{2}^{3})$, and the $k$ eigenvectors in
$\mathcal{O}(kN)$ operations. For $N_{1}\approx N_{2}\approx\sqrt{N}$, exact
sampling thus only costs $\mathcal{O}(N^{3/2}+Nk^{3})$. If $L=L_{1}\otimes
L_{2}\otimes L_{3}$, the same reasoning shows that exact sampling becomes
linear in $N$, only requiring $\mathcal{O}(Nk^{3})$ operations.

One can also resort to MCMC sampling; for instance such a sampler was
considered in [13] (though with an incorrect mixing time analysis). The
results of [20] hold only for $k$-DPPs, but suggest their MCMC sampler may
possibly take $O(N^{2}\log(N/\epsilon))$ time for full DPPs, which is
impractical. Nevertheless if one develops faster MCMC samplers, they should
also be able to profit from the Kronecker product structure offered by
KronDpp.

##  5 Experimental results

In order to validate our learning algorithm, we compared KrK-Picard to Joint-
Picard and to the Picard iteration (Picard) on multiple real and synthetic
datasets.111All experiments were repeated 5 times and averaged, using MATLAB
on a Linux Mint system with 16GB of RAM and an i7-4710HQ CPU @ 2.50GHz.

###  5.1 Synthetic tests

All three algorithms were used to learn from synthetic data drawn from a
“true” kernel. The sub-kernels were initialized by $L_{i}=X^{\top}X$, with
$X$’s coefficients drawn uniformly from $[0,\sqrt{2}]$; for Picard, $L$ was
initialized with $L_{1}\otimes L_{2}$.

(a) $N_{1}=N_{2}=50$

(b) $N_{1}=N_{2}=100$

(c) $N_{1}=100,N_{2}=500$

Figure 1: $a=1$; the thin dotted lines indicated the standard deviation from
the mean.

For Figures 1(a) and 1(b), training data was generated by sampling 100 subsets
from the true kernel with sizes uniformly distributed between 10 and 190.

To evaluate KrK-Picard on matrices too large to fit in memory and with large
$\kappa$, we drew samples from a $50\cdot 10^{3}\\!\times\\!50\cdot 10^{3}$
kernel of rank $1,000$ (on average $|Y_{i}|\approx 1,000$), and learned the
kernel stochastically (only KrK-Picard could be run due to the memory
requirements of other methods); the likelihood drastically improves in only
two steps (Fig.1(c)).

As shown in Figures 1(a) and 1(b), KrK-Picard converges significantly faster
than Picard, especially for large values of $N$. However, although Joint-
Picard also increases the log-likelihood at each iteration, it converges much
slower and has a high standard deviation, whereas the standard deviations for
Picard and KrK-Picard are barely noticeable. For these reasons, we drop the
comparison to Joint-Picard in the subsequent experiments.

###  5.2 Small-scale real data: baby registries

We compared KrK-Picard to Picard and EM [10] on the baby registry dataset
(described in-depth in [10]), which has also been used to evaluate other DPP
learning algorithms [10, 25, 9]. The dataset contains 17 categories of baby-
related products obtained from Amazon. We learned kernels for the 6 largest
categories ($N=100$); in this case, Picard is sufficiently efficient to be
prefered to KrK-Picard; this comparison serves only to evaluate the quality of
the final kernel estimates.

The initial marginal kernel $K$ for EM was sampled from a Wishart distribution
with $N$ degrees of freedom and an identity covariance matrix, then scaled by
$1/N$; for Picard, $L$ was set to $K(I-K)^{-1}$; for KrK-Picard, $L_{1}$ and
$L_{2}$ were chosen (as in Joint-Picard) by minimizing $\|L-L_{1}\otimes
L_{2}\|$. Convergence was determined when the objective change dipped below a
threshold $\delta$. As one EM iteration takes longer than one Picard iteration
but increases the likelihood more, we set
$\delta_{\textsc{Pic}}=\delta_{\textsc{KrK}}=10^{-4}$ and
$\delta_{\textsc{EM}}=10^{-5}$.

The final log-likelihoods are shown in Table 1(b); we set the step-sizes to
their largest possible values, i.e. $a_{\textsc{Pic}}=1.3$ and
$a_{\textsc{KrK}}=1.8$. Table 1(b) shows that KrK-Picard obtains comparable,
albeit slightly worse log-likelihoods than Picard and EM, which confirms that
for tractable $N$, the better modeling capability of full kernels make them
preferable to KronDpps.

Table 1: Final log-likelihoods for each large category of the baby registries
dataset

(a) Training set Category | EM | Picard | KrK-Picard
---|---|---|---
apparel | -10.1 | -10.2 | -10.7
bath | -8.6 | -8.8 | -9.1
bedding | -8.7 | -8.8 | -9.3
diaper | -10.5 | -10.7 | -11.1
feeding | -12.1 | -12.1 | -12.5
gear | -9.3 | -9.3 | -9.6

(b) Test set Category | EM | Picard | KrK-Picard
---|---|---|---
apparel | -10.1 | -10.2 | -10.7
bath | -8.6 | -8.8 | -9.1
bedding | -8.7 | -8.8 | -9.3
diaper | -10.6 | -10.7 | -11.2
feeding | -12.2 | -12.2 | -12.6
gear | -9.2 | -9.2 | -9.5

###  5.3 Large-scale real dataset: GENES

Finally, to evaluate KrK-Picard on large matrices of real-world data, we train
it on data from the GENES [3] dataset (which has also been used to evaluate
DPPs in [19, 3]). This dataset consists in 10,000 genes, each represented by
331 features corresponding to the distance of a gene to hubs in the BioGRID
gene interaction network.

Table 2: Average runtime and performance on the GENES dataset for $N_{1}=N_{2}=100$ | Picard | KrK-Picard |  KrK-Picard (stochastic)
---|---|---|---
Average runtime |  161.5 $\pm$ 17.7 s |  8.9 $\pm$ 0.2 s |  1.2 $\pm$ 0.02 s
NLL increase (1st iter.) | $(2.81\pm 0.03)\cdot 10^{4}$ | $(2.96\pm 0.02)\cdot 10^{4}$ | $(3.13\pm 0.04)\cdot 10^{4}$

We construct a ground truth Gaussian DPP kernel on the GENES dataset and use
it to obtain 100 training samples with sizes uniformly distributed between 50
and 200 items. Similarly to the synthetic experiments, we initialized KrK-
Picard’s kernel by setting $L_{i}=X_{i}^{\top}X_{i}$ where $X_{i}$ is a random
matrix of size $N_{1}\times N_{1}$; for Picard, we set the initial kernel
$L=L_{1}\otimes L_{2}$.

(a) Non-stochastic learning

(b) Stochastic vs. non-stochastic

Figure 2: $n=150$, $a=1$.

Figure 2 shows the performance of both algorithms. As with the synthetic
experiments, KrK-Picard converges much faster; stochastic updates increase its
performance even more, as shown in Fig. 2(b). Average runtimes and speed-up
are given in Table 2: KrK-Picard runs almost an order of magnitude faster than
Picard, and stochastic updates are more than two orders of magnitude faster,
while providing slightly larger initial increases to the log-likelihood.

##  6 Conclusion and future work

We introduced KronDpps, a variant of DPPs with kernels structured as the
Kronecker product of $m$ smaller matrices, and showed that typical operations
over DPPs such as sampling and learning the kernel from data can be made
efficient for KronDpps on previously untractable ground set sizes.

By carefully leveraging the properties of the Kronecker product, we derived
for $m=2$ a low-complexity algorithm to learn the kernel from data which
guarantees positive iterates and a monotonic increase of the log-likelihood,
and runs in $\mathcal{O}(n\kappa^{3}+N^{2})$ time. This algorithm provides
even more significant speed-ups and memory gains in the stochastic case,
requiring only $\mathcal{O}(N^{3/2}+N\kappa^{2})$ time and
$\mathcal{O}(N+\kappa^{2})$ space. Experiments on synthetic and real data
showed that KronDpps can be learned efficiently on sets large enough that $L$
does not fit in memory.

While discussing learning the kernel, we showed that $L_{1}$ and $L_{2}$
cannot be updated simultaneously in a CCCP-style iteration since $g$ is not
convex over $(S_{1},S_{2})$. However, it can be shown that $g$ is geodesically
convex over the Riemannian manifold of positive definite matrices, which
suggests that deriving an iteration which would take advantage of the
intrinsic geometry of the problem may be a viable line of future work.

KronDpps also enable fast sampling, in $\mathcal{O}(N^{3/2}+Nk^{3})$
operations when using two sub-kernels and in $\mathcal{O}(Nk^{3})$ when using
three sub-kernels; this allows for exact sampling at comparable or even better
costs than previous algorithms for approximate sampling. However, as we
improve computational efficiency $L$, the subset size $k$ becomes limiting,
due to the $\mathcal{O}(Nk^{3})$ cost of sampling and learning. A necessary
line of future work to allow for truly scalable DPPs is thus to overcome this
computational bottleneck.

## References

  * Affandi et al. [2013] R. Affandi, A. Kulesza, E. Fox, and B. Taskar.  Nyström approximation for large-scale Determinantal Point Processes.  In _Artificial Intelligence and Statistics (AISTATS)_ , 2013.
  * Affandi et al. [2014] R. Affandi, E. Fox, R. Adams, and B. Taskar.  Learning the parameters of Determinantal Point Process kernels.  In _ICML_ , 2014.
  * Batmanghelich et al. [2014] N. K. Batmanghelich, G. Quon, A. Kulesza, M. Kellis, P. Golland, and L. Bornn.  Diversifying sparsity using variational determinantal point processes.  _arXiv:1411.6307_ , 2014.
  * Bhatia [2007] R. Bhatia.  _Positive Definite Matrices_.  Princeton University Press, 2007.
  * Borodin [2009] A. Borodin.  Determinantal point processes.  _arXiv:0911.1153_ , 2009.
  * Chao et al. [2015] W. Chao, B. Gong, K. Grauman, and F. Sha.  Large-margin determinantal point processes.  In _Uncertainty in Artificial Intelligence (UAI)_ , 2015.
  * Decreusefond et al. [2015] L. Decreusefond, I. Flint, N. Privault, and G. L. Torrisi.  Determinantal point processes, 2015.
  * Flaxman et al. [2015] S. Flaxman, A. Wilson, D. Neill, H. Nickisch, and A. Smola.  Fast Kronecker inference in Gaussian processes with non-Gaussian likelihoods.  In _ICML_ , pages 607–616, 2015.
  * Gartrell et al. [2016] M. Gartrell, U. Paquet, and N. Koenigstein.  Low-rank factorization of determinantal point processes for recommendation.  _arXiv:1602.05436_ , 2016.
  * Gillenwater et al. [2014] J. Gillenwater, A. Kulesza, E. Fox, and B. Taskar.  Expectation-Maximization for learning Determinantal Point Processes.  In _NIPS_ , 2014.
  * Goldschmidt et al. [1994] O. Goldschmidt, D. Nehme, and G. Yu.  Note: On the set-union knapsack problem.  _Naval Research Logistics_ , 41:833–842, 1994.
  * Hough et al. [2006] J. B. Hough, M. Krishnapur, Y. Peres, and B. Virág.  Determinantal processes and independence.  _Probability Surveys_ , 3(206–229):9, 2006.
  * Kang [2013] B. Kang.  Fast determinantal point process sampling with application to clustering.  In _Advances in Neural Information Processing Systems 26_ , pages 2319–2327. Curran Associates, Inc., 2013.
  * Krause et al. [2008] A. Krause, A. Singh, and C. Guestrin.  Near-optimal sensor placements in Gaussian processes: theory, efficient algorithms and empirical studies.  _JMLR_ , 9:235–284, 2008.
  * Kulesza [2013] A. Kulesza.  _Learning with Determinantal Point Processes_.  PhD thesis, University of Pennsylvania, 2013.
  * Kulesza and Taskar [2011] A. Kulesza and B. Taskar.  k-DPPs: Fixed-size Determinantal Point Processes.  In _ICML_ , 2011.
  * Kulesza and Taskar [2012] A. Kulesza and B. Taskar.  _Determinantal Point Processes for machine learning_ , volume 5.  Foundations and Trends in Machine Learning, 2012.
  * Lavancier et al. [2015] F. Lavancier, J. Møller, and E. Rubak.  Determinantal point process models and statistical inference.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 77(4):853–877, 2015.
  * Li et al. [2015] C. Li, S. Jegelka, and S. Sra.  Efficient sampling for k-determinantal point processes.  _arXiv:1509.01618_ , 2015.
  * Li et al. [2016] C. Li, S. Jegelka, and S. Sra.  Fast DPP sampling for Nyström with application to kernel methods.  _arXiv:1603.06052_ , 2016.
  * Lin and Bilmes [2012] H. Lin and J. Bilmes.  Learning mixtures of submodular shells with application to document summarization.  In _Uncertainty in Artificial Intelligence (UAI)_ , 2012.
  * Loan and Pitsianis [1993] C. V. Loan and N. Pitsianis.  Approximation with kronecker products.  In _Linear Algebra for Large Scale and Real Time Applications_ , pages 293–314. Kluwer Publications, 1993.
  * Lyons [2003] R. Lyons.  Determinantal probability measures.  _Publications Mathématiques de l’Institut des Hautes Études Scientifiques_ , 98(1):167–212, 2003.
  * Macchi [1975] O. Macchi.  The coincidence approach to stochastic point processes.  _Adv. Appl. Prob._ , 7(1), 1975.
  * Mariet and Sra [2015] Z. Mariet and S. Sra.  Fixed-point algorithms for learning determinantal point processes.  In _ICML_ , 2015.
  * Mariet and Sra [2016] Z. Mariet and S. Sra.  Diversity networks.  _Int. Conf. on Learning Representations (ICLR)_ , 2016.  URL arXiv:1511.05077.
  * Martens and Grosse [2015] J. Martens and R. B. Grosse.  Optimizing neural networks with Kronecker-factored approximate curvature.  In _ICML_ , 2015.
  * Wu et al. [2005] G. Wu, Z. Zhang, and E. Y. Chang.  Kronecker factorization for speeding up kernel machines.  In _SIAM Data Mining (SDM)_ , pages 611–615, 2005.
  * Yuille and Rangarajan [2002] A. L. Yuille and A. Rangarajan.  The concave-convex procedure (cccp).  In _Advances in Neural Information Processing Systems 14_ , pages 1033–1040. MIT Press, 2002.
  * Zhang et al. [2015] X. Zhang, F. X. Yu, R. Guo, S. Kumar, S. Wang, and S.-F. Chang.  Fast orthogonal projection based on kronecker product.  In _ICCV_ , 2015.
  * Zhou et al. [2010] T. Zhou, Z. Kuscsik, J.-G. Liu, M. Medo, J. R. Wakeling, and Y.-C. Zhang.  Solving the apparent diversity-accuracy dilemma of recommender systems.  _PNAS_ , 107(10):4511–4515, 2010.

Appendix: Kronecker Determinantal Point Processes

##  Appendix A Proof of Prop. 3.1

We use ‘$\operatorname{vec}$’ to denote the operator that stacks columns of a
matrix to form a vector; conversely, ‘$\operatorname{mat}$’ takes a vector
with $k^{2}$ coefficients and returns a $k\times k$ matrix.

Let $L=L_{1}\otimes L_{2}$, $S_{1}=L_{1}^{-1},S_{2}=L_{2}^{-1}$ and
$S=S_{1}\otimes S_{2}=L^{-1}$. We note $E_{ij}$ the matrix with all zeros
except for a 1 at position $(i,j)$, its size being clear from context. We wish
to solve

$\nabla\\!f_{2}(X)=-\nabla\\!g_{2}(S_{1})\quad\text{ and
}\nabla\\!f_{1}(X)=-\nabla\\!g_{1}(S_{2})$ (10)

It follows from the fact that

$\log\det(S_{1}\otimes S_{2})=N_{2}\log\det S_{1}+N_{1}\log\det S_{2}$

that $\nabla\\!f_{S_{2}}(X)=N_{2}X^{-1}$ and
$\nabla\\!f_{S_{1}}(X)=N_{1}X^{-1}$. Moreover, we know that

$\displaystyle\nabla\\!g(S)$
$\displaystyle=-(I+S)^{-1}-S^{-1}\frac{1}{n}\sum\nolimits_{i}U_{i}(U_{i}^{\top}S^{-1}U_{i})^{-1}U_{i}S^{-1}$
$\displaystyle=-S^{-1}-S^{-1}\left(\frac{1}{n}\sum\nolimits_{i}U_{i}(U_{i}^{\top}S^{-1}U_{i})^{-1}U_{i}-(I+S^{-1})^{-1}\right)S^{-1}$
$\displaystyle=-(L+L\Delta L).$

The Jacobian of $S_{1}\rightarrow S_{1}\otimes S_{2}$ is given by
$J=\begin{pmatrix}\operatorname{vec}(E_{11}\otimes
S_{2}),\ldots,\operatorname{vec}(E_{N_{1}N_{1}}\otimes S_{2})\end{pmatrix}$.
Hence,

$\displaystyle\nabla\\!f_{1}(X)_{ij}=-(\nabla\\!g_{1}(S_{1}))_{ij}$
$\displaystyle\iff N_{2}X^{-1}_{ij}=(J^{\top}\operatorname{vec}(-\nabla
g(S)))_{ij}$ $\displaystyle\iff
N_{2}X^{-1}_{ij}=\operatorname{vec}(E_{ij}\otimes
S_{2})^{\top}\operatorname{vec}(L+L\Delta L)$ $\displaystyle\iff
N_{2}X^{-1}_{ij}=\operatorname{Tr}((E_{ij}\otimes S_{2})(L+L\Delta L))$
$\displaystyle\iff N_{2}X^{-1}_{ij}=\operatorname{Tr}(S_{2}(L+L\Delta
L)_{(ij)})$ $\displaystyle\iff
N_{2}X^{-1}_{ij}=\operatorname{Tr}\left(((I\otimes S_{2})(L+L\Delta
L))_{(ij)}\right)$

The last equivalence is simply the result of indices manipulation. Thus, we
have

$\nabla f_{2}(X)=-\nabla g_{2}(S_{1})\iff
X^{-1}=\frac{1}{N_{2}}\operatorname{Tr}_{1}\\!\left((I\otimes S_{2})(L+L\Delta
L)\right)$

Similarly, by setting $J^{\prime}=\left(\operatorname{vec}(S_{1}\otimes
E_{11}),\ldots,\operatorname{vec}(S_{1}\otimes E_{N_{1}N_{1}})\right)$, we
have that

$\displaystyle\nabla\\!f_{2}(X)_{ij}=-(\nabla\\!g_{2}(S_{2}))_{ij}$
$\displaystyle\iff N_{1}X^{-1}_{ij}=(J^{\prime\top}\operatorname{vec}(-\nabla
g(S)))_{ij}$ $\displaystyle\iff
N_{1}X^{-1}_{ij}=\operatorname{vec}(S_{1}\otimes
E_{ij})^{\top}\operatorname{vec}(L+L\Delta L)$ $\displaystyle\iff
N_{1}X^{-1}_{ij}=\operatorname{Tr}((S_{1}\otimes E_{ij})(L+L\Delta L))$
$\displaystyle\iff
N_{1}X^{-1}_{ij}=\left(\sum\nolimits_{k,\ell=1}^{N_{1}}{S_{1}}_{k\ell}(L+L\Delta
L)_{(\ell k)}\right)_{ij}$ $\displaystyle\iff
N_{1}X^{-1}_{ij}=\left(\sum\nolimits_{\ell=1}^{N_{1}}((S_{1}\otimes
I)(L+L\Delta L))_{(\ell\ell)}\right)_{ij}$

Hence,

$\nabla f_{S_{1}}(X)=-\nabla g_{S_{1}}(S_{2})\iff
X^{-1}=\frac{1}{N_{1}}\operatorname{Tr}_{2}\left((S_{1}\otimes I)(L+L\Delta
L)\right),$

which proves Prop. 3.1. ∎

##  Appendix B Efficient updates for KrK-Picard

The updates to $L_{1}$ and $L_{2}$ are obtained efficiently through different
methods; hence, the proof to Thm. 3.3 is split into two sections. We write

$\Theta=\frac{1}{n}\sum_{i=1}^{n}U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}\quad\text{
(or }\Theta=U_{i}L_{Y_{i}}^{-1}U_{i}^{\top}\text{ for stochastic updates)}$

so that $\Delta=\Theta-(I+L)^{-1}$. Recall that $(A\otimes B)_{(ij)}=a_{ij}B$.

###  B.1 Updating $L_{1}$

We wish to compute $X=\operatorname{Tr}_{1}\\!\left((I\otimes
L_{2}^{-1})(L\Delta L)\right)$ efficiently. We have

$\displaystyle X_{ij}$ $\displaystyle=\operatorname{Tr}\left[((I\otimes
L_{2}^{-1})(L\Delta L))_{(ij)}\right]$
$\displaystyle=\operatorname{Tr}\left[L_{2}^{-1}(L\Delta L)_{(ij)}\right]$
$\displaystyle=\operatorname{Tr}\left[L_{2}^{-1}\sum\nolimits_{k,\ell=1}^{N_{1}}L_{(ik)}\Delta_{(k\ell)}L_{(\ell
j)}\right]$
$\displaystyle=\sum\nolimits_{k,\ell=1}^{N_{1}}{L_{1}}_{ik}{L_{1}}_{\ell
j}\operatorname{Tr}(L_{2}^{-1}L_{2}\Delta_{(k\ell)}L_{2})$
$\displaystyle=\sum\nolimits_{k,\ell=1}^{N_{1}}{L_{1}}_{ik}{L_{1}}_{\ell
j}\underbrace{\operatorname{Tr}(\Theta_{(k\ell)}L_{2})}_{A_{k\ell}}-\underbrace{\operatorname{Tr}((I+L)^{-1}_{(k\ell)}L_{2})}_{B_{k\ell}}$
$\displaystyle=(L_{1}AL_{1}-L_{1}BL_{1})_{ij}.$

The $N_{1}\times N_{1}$ matrix $A$ can be computed in
$\mathcal{O}(n\kappa^{3}+N_{1}^{2}N_{2}^{2})$ simply by pre-computing $\Theta$
in $\mathcal{O}(n\kappa^{3})$ and then computing all $N_{1}^{2}$ traces in
$\mathcal{O}(N_{2}^{2})$ time. When doing stochastic updates for which
$\Theta$ is sparse with only $\kappa^{2}$ non-zero coefficients, computing $A$
can be done in $\mathcal{O}(N_{1}^{2}\kappa^{2}+\kappa^{3})$.

By diagonalizing $L_{1}=P_{1}D_{1}P_{1}^{\top}$ and
$L_{2}=P_{2}D_{2}P_{2}^{\top}$, we have $(I+L)^{-1}=PDP^{\top}$ with
$P=P_{1}\otimes P_{2}$ and $D=(I+D_{1}\otimes D_{2})^{-1}$.
$P_{1},P_{2},D_{1},D_{2}$ and $D$ can all be obtained in
$\mathcal{O}(N_{1}^{3}+N_{2}^{3}+N_{1}N_{2})$ as a consequence of Prop. 2.1.
Then

$\displaystyle B_{ij}$
$\displaystyle=\operatorname{Tr}((I+L)^{-1}_{(ij)}L_{2})$
$\displaystyle=\sum_{k}\operatorname{Tr}(P_{(ik)}D_{(kk)}P^{\top}_{(kj)}L_{2})$
$\displaystyle=\sum_{k}{P_{1}}_{ik}{P_{1}}_{jk}\operatorname{Tr}(P_{2}D_{(kk)}P_{2}^{\top}P_{2}D_{2}P_{2}^{\top})$
$\displaystyle=\sum_{k}{P_{1}}_{ik}{P_{1}}_{jk}\underbrace{\operatorname{Tr}(D_{(kk)}D_{2})}_{\alpha_{k}}.$

Let $\widehat{D}=\text{diag}(\alpha_{1},\ldots,\alpha_{N_{1}})$, which can be
computed in $\mathcal{O}(N_{1}N_{2})$. Then
$L_{1}BL_{1}=P_{1}D_{1}\widehat{D}D_{1}P_{1}$ is computable in
$\mathcal{O}(N_{1}^{3}+N_{2}^{3})$.

Overall, the update to $L_{1}$ can be computed in
$\mathcal{O}(n\kappa^{3}+N_{1}^{2}N_{2}^{2}+N_{1}^{3}+N_{2}^{3})$, or in
$\mathcal{O}(N_{1}^{2}\kappa^{2}+\kappa^{3}+N_{1}^{3}+N_{2}^{3})$ if the
updates are stochastic. Moreover, if $\Theta$ is sparse with only $z$ non-zero
coefficients (for stochastic updates $z=\kappa$), $A$ can be computed in
$\mathcal{O}(\kappa^{2})$ space, leading to an overall
$\mathcal{O}(z^{2}+N_{1}^{2}+N_{2}^{2})$ memory cost.

###  B.2 Updating $L_{2}$

We wish to compute $X=\operatorname{Tr}_{2}\\!\left[(L_{1}^{-1}\otimes
I)(L\Delta L)\right]$ efficiently.

$\displaystyle X$
$\displaystyle=\sum\nolimits_{i=1}^{N_{1}}\left((L_{1}^{-1}\otimes I)(L\Delta
L)\right)_{(ii)}$ $\displaystyle=\sum\nolimits_{i=1}^{N_{1}}\left((I\otimes
L_{2})(\Theta-(I+L)^{-1})(L_{1}\otimes L_{2})\right)_{(ii)}$
$\displaystyle=\sum\nolimits_{i,j=1}^{N_{1}}{L_{1}}_{ij}L_{2}\Theta_{(ij)}L_{2}-\sum_{i=1}^{N_{1}}((I\otimes
L_{2})(I+L)^{-1}(L_{1}\otimes L_{2}))_{(ii)}$
$\displaystyle=\underbrace{L_{2}\sum\nolimits_{i,j=1}^{N_{1}}{L_{1}}_{ij}\Theta_{(ij)}L_{2}}_{A}-\underbrace{\sum\nolimits_{i=1}^{N_{1}}((I\otimes
L_{2})(I+L)^{-1}(L_{1}\otimes L_{2}))_{(ii)}}_{B}$

$A$ can be computed in
$\mathcal{O}(n\kappa^{3}+N_{1}^{2}N_{2}^{2}+N_{2}^{3})$. As before, when doing
stochastic updates $A$ can be computed in
$\mathcal{O}(N_{1}^{2}\kappa^{2}+\kappa^{3}+N_{2}^{3})$ time and
$\mathcal{O}(N_{2}^{2}+N_{1}^{2}+\kappa^{2})$ space due to the sparsity of
$\Theta$.

Regarding $B$, as all matrices commute, we can write

$(I\otimes L_{2})(I+L)^{-1}(L_{1}\otimes L_{2})=(P_{1}\otimes
P_{2})\Lambda(P_{1}\otimes P_{2})$

where $\Lambda=(I\otimes D_{2})(I+D_{1}\otimes D_{2})^{-1}(D_{1}\otimes
D_{2})$ is diagonal and is obtained in
$\mathcal{O}(N_{1}^{3}+N_{2}^{3}+N_{1}N_{2})$. Moreover,

$B=\sum\nolimits_{i=1}^{N_{1}}(P\Lambda
P^{\top})_{(ii)}=P_{2}\left(\sum\nolimits_{i,k=1}^{N_{1}}{P_{1}}_{ik}\Lambda_{(kk)}{P_{1}}_{ik}\right)P_{2}^{\top},$

which allows us to compute $B$ in
$\mathcal{O}(N_{1}^{2}N_{2}+N_{2}^{3}+N_{1}^{3})$ total.

Overall, we can obtain $X$ in
$\mathcal{O}(n\kappa^{3}+N_{1}^{2}N_{2}^{2}+N_{1}^{3}+N_{2}^{3})$ or in
$\mathcal{O}(N_{1}^{2}\kappa^{2}+N_{1}^{2}N_{2}+N_{1}^{3}+N_{2}^{3})$ for
stochastic updates, in which case only
$\mathcal{O}(N_{1}^{2}+N_{2}^{2}+\kappa^{2})$ space is necessary.

##  Appendix C Proof of validity for joint updates

In order to minimize the number of matrix multiplications, we equivalently
(due to the properties of the Frobenius norm) minimize the equation

$\|L^{-1}+\Delta-X\otimes Y\|^{2}_{F}$ (11)

and set $\begin{cases}L_{1}^{\prime}\leftarrow L_{1}XL_{1}\\\
L_{2}^{\prime}\leftarrow L_{2}YL_{2}.\end{cases}$.

######  Theorem C.1.

Let $L\succ 0$. Define
$R:=[\operatorname{vec}(L_{(11)})^{\top};\ldots;\operatorname{vec}(L_{(N_{1}N_{1})})^{\top}]_{i,j=1}^{N_{1}}\in\mathbb{R}^{N_{1}N_{1}\times
N_{2}N_{2}}$.

Suppose that $R$ has an eigengap between its largest singular value and the
next, and let $u,v,\sigma$ be the first singular vectors and value of $R$. Let
$U=\operatorname{mat}(u)$ and $V=\operatorname{mat}(v)$. Then $U$ and $V$ are
either both positive definite or negative definite.

Moreover, for any value $\alpha\neq 0$, the pair $(\alpha U,\sigma/\alpha V)$
minimizes $\|L-X\otimes Y\|^{2}_{F}$.

The proof is a consequence of [22, Thm. 11]. This shows that if $L$ is
initially positive definite, setting the sign of $\alpha$ based on whether $U$
and $V$ are positive or negative definite222This can easily be done simply by
checking the sign of the first diagonal coefficient of $U$, which will be
positive if and only if $U\succ 0$., and updating

$\begin{cases}L_{1}\leftarrow\alpha\,L_{1}UL_{1}\\\
L_{2}\leftarrow\sigma/\alpha\,L_{2}VL_{2}\end{cases}$

maintains positive definite iterates. Given that if $L_{1}\succ 0$ and
$L_{2}\succ 0$, $L_{1}\otimes L_{2}\succ 0$, a simple induction then shows
that by choosing an initial kernel estimate $L\succ 0$, subsequent values of
$L$ will remain positive definite.

By choosing $\alpha$ such that the new estimates $L_{1}$ and $L_{2}$ verify
$\|L_{1}\|=\|L_{2}\|$, we verify all the conditions of Eq. 8.

###  C.1 Algorithm for joint updates

Theorem C.1 leads to a straightforward iteration for learning matrices $L_{1}$
and $L_{2}$ based on the decomposition of the Picard estimate as a Kronecker
product.

Algorithm 3 Joint-Picard iteration

Input: Matrices $L_{1},L_{2}$, training set $T$, step-size $a\geq 1$.

for $i=1$ to maxIter do

$U,\sigma,V\leftarrow$ power_method$(L^{-1}+\Delta)$ to obtain the first
singular value and vectors of matrix $R$.

$\alpha\leftarrow\text{sgn}(U_{11})\sqrt{\sigma\|L_{2}VL_{2}\|/\|L_{1}UL_{1}\|}$

$L_{1}\leftarrow L_{1}+a(\alpha\,L_{1}UL_{1}-L_{1})$

$L_{2}\leftarrow L_{2}+a(\sigma/\alpha\,L_{2}VL_{2})$

end forreturn $(L_{1},L_{2})$