Expectation-Maximization for Learning Determinantal Point Processes

Introduction

Subset selection is a core task in many real-world applications. For example, in product recommendation we typically want to choose a small set of products from a large collection; many other examples of subset selection tasks turn up in domains like document summarization Lin and Bilmes (2012); Kulesza and Taskar (2012), sensor placement Krause et al. (2008); Krause and Guestrin (2005), image search Kulesza and Taskar (2011b); Affandi et al. (2014), and auction revenue maximization Dughmi et al. (2009), to name a few. In these applications, a good subset is often one whose individual items are all high-quality, but also all distinct. For instance, recommended products should be popular, but they should also be diverse to increase the chance that a user finds at least one of them interesting. Determinantal point processes (DPPs) offer one way to model this tradeoff; a DPP defines a distribution over all possible subsets of a ground set, and the mass it assigns to any given set is a balanced measure of that set’s quality and diversity.

Originally discovered as models of fermions Macchi (1975), DPPs have recently been effectively adapted for a variety of machine learning tasks Affandi et al. (2014); Snoek et al. (2013); Kang (2013); Affandi et al. (2013a); Shah and Ghahramani (2013); Affandi et al. (2013b); Gillenwater et al. (2012a); Zou and Adams (2013); Affandi et al. (2012); Gillenwater et al. (2012b); Kulesza and Taskar (2011a, b, 2010). They offer attractive computational properties, including exact and efficient normalization, marginalization, conditioning, and sampling Hough et al. (2006). These properties arise in part from the fact that a DPP can be compactly parameterized by an $N\times N$ positive semi-definite matrix $L$. Unfortunately, though, learning $L$ from example subsets by maximizing likelihood is conjectured to be NP-hard (Kulesza, 2012, Conjecture 4.1). While gradient ascent can be applied in an attempt to approximately optimize the likelihood objective, we show later that it requires a projection step that often produces degenerate results.

For this reason, in most previous work only partial learning of $L$ has been attempted. Kulesza and Taskar (2011a) showed that the problem of learning a scalar weight for each row of $L$ is a convex optimization problem. This amounts to learning what makes an item high-quality, but does not address the issue of what makes two items similar. Kulesza and Taskar (2011b) explored a different direction, learning weights for a linear combination of DPPs with fixed $L$s. This works well in a limited setting, but requires storing a potentially large set of kernel matrices, and the final distribution is no longer a DPP, which means that many attractive computational properties are lost. Affandi et al. (2014) proposed as an alternative that one first assume $L$ takes on a particular parametric form, and then sample from the posterior distribution over kernel parameters using Bayesian methods. This overcomes some of the disadvantages of Kulesza and Taskar (2011b)’s $L$-ensemble method, but does not allow for learning an unconstrained, non-parametric $L$.

The learning method we propose in this paper differs from those of prior work in that it does not assume fixed values or restrictive parameterizations for $L$, and exploits the eigendecomposition of $L$. Many properties of a DPP can be simply characterized in terms of the eigenvalues and eigenvectors of $L$, and working with this decomposition allows us to develop an expectation-maximization (EM) style optimization algorithm. This algorithm negates the need for the problematic projection step that is required for naive gradient ascent to maintain positive semi-definiteness of $L$. As the experiments show, a projection step can sometimes lead to learning a nearly diagonal $L$, which fails to model the negative interactions between items. These interactions are vital, as they lead to the diversity-seeking nature of a DPP. The proposed EM algorithm overcomes this failing, making it more robust to initialization and dataset changes. It is also asymptotically faster than gradient ascent.

Background

An alternative representation of a DPP is given by the marginal kernel: $K=L(L+I)^{-1}$. The $L$-$K$ relationship can also be written in terms of their eigendecompositons. $L$ and $K$ share the same eigenvectors $\boldsymbol{v}$, and an eigenvalue $\lambda_{i}$ of $K$ corresponds to an eigenvalue $\lambda_{i}/(1-\lambda_{i})$ of $L$:

Clearly, if $L$ if PSD then $K$ is as well, and the above equations also imply that the eigenvalues of $K$ are further restricted to be $\leq 1$. $K$ is called the marginal kernel because, for any set $Y\sim\mathcal{P}$ and for every $A\subseteq\mathcal{Y}$:

We can also write the exact (non-marginal, normalized) probability of a set $Y\sim\mathcal{P}$ in terms of $K$:

where $I_{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{Y}}$ is the identity matrix with entry $(i,i)$ zeroed for items $i\in Y$ (Kulesza, 2012, Equation 3.69). In what follows we use the $K$-based formula for $\mathcal{P}(Y)$ and learn the marginal kernel $K$. This is equivalent to learning $L$, as Equation (1) can be applied to convert from $K$ to $L$.

Learning algorithms

In our learning setting the input consists of $n$ example subsets, $\{Y_{1},\ldots,Y_{n}\}$, where $Y_{i}\subseteq\{1,\ldots,N\}$ for all $i$. Our goal is to maximize the likelihood of these example sets. We first describe in Section 3.1 a naive optimization procedure: projected gradient ascent on the entries of the marginal matrix $K$, which will serve as a baseline in our experiments. We then develop an EM method: Section 3.2 changes variables from kernel entries to eigenvalues and eigenvectors (introducing a hidden variable in the process), Section 3.3 applies Jensen’s inequality to lower-bound the objective, and Sections 3.4 and 3.5 outline a coordinate ascent procedure on this lower bound.

The log-likelihood maximization problem, based on Equation (3), is:

where the first constraint ensures that $K$ is PSD and the second puts an upper limit of $1$ on its eigenvalues. Let $\mathcal{L}(K)$ represent this log-likelihood objective. Its partial derivative with respect to $K$ is easy to compute by applying a standard matrix derivative rule (Petersen and Pedersen, 2012, Equation 57):

Thus, projected gradient ascent Levitin and Polyak (1966) is a viable, simple optimization technique. Algorithm 1 outlines this method, which we refer to as K-Ascent (KA). The initial $K$ supplied as input to the algorithm can be any PSD matrix with eigenvalues $\leq 1$. The first part of the projection step, $\max(\boldsymbol{\lambda},0)$, chooses the closest (in Frobenius norm) PSD matrix to $Q$ (Henrion and Malick, 2011, Equation 1). The second part, $\min(\boldsymbol{\lambda},1)$, caps the eigenvalues at $1$. (Notice that only the eigenvalues have to be projected; $K$ remains symmetric after the gradient step, so its eigenvectors are already guaranteed to be real.)

Unfortunately, the projection can take us to a poor local optima. To see this, consider the case where the starting kernel $K$ is a poor fit to the data. In this case, a large initial step size $\eta$ will probably be accepted; even though such a step will likely result in the truncation of many eigenvalues at , the resulting matrix will still be an improvement over the poor initial $K$. However, with many zero eigenvalues, the new $K$ will be near-diagonal, and, unfortunately, Equation (5) dictates that if the current $K$ is diagonal, then its gradient is as well. Thus, the KA algorithm cannot easily move to any highly non-diagonal matrix. It is possible that employing more complex step-size selection mechanisms could alleviate this problem, but the EM algorithm we develop in the next section will negate the need for these entirely.

The EM algorithm we develop also has an advantage in terms of asymptotic runtime. The computational complexity of KA is dominated by the matrix inverses of the $\mathcal{L}$ derivative, each of which requires $O(N^{3})$ operations, and by the eigendecomposition needed for the projection, also $O(N^{3})$. The overall runtime of KA, assuming $T_{1}$ iterations until convergence and an average of $T_{2}$ iterations to find a step size, is $O(T_{1}nN^{3}+T_{1}T_{2}N^{3})$. As we will show in the following sections, the overall runtime of the EM algorithm is $O(T_{1}nNk^{2}+T_{1}T_{2}N^{3})$, which can be substantially better than KA’s runtime for $k\ll N$.

2 Eigendecomposing

Eigendecomposition is key to many core DPP algorithms such as sampling and marginalization. This is because the eigendecomposition provides an alternative view of the DPP as a generative process, which often leads to more efficient algorithms. Specifically, sampling a set $Y$ can be broken down into a two-step process, the first of which involves generating a hidden variable $J\subseteq\{1,\ldots,N\}$ that codes for a particular set of $K$’s eigenvectors. We review this process below, then exploit it to develop an EM optimization scheme.

Suppose $K=V\Lambda V^{\top}$ is an eigendecomposition of $K$. Let $V^{J}$ denote the submatrix of $V$ containing only the columns corresponding to the indices in a set $J\subseteq\{1,\ldots,N\}$. Consider the corresponding marginal kernel, with all selected eigenvalues set to $1$:

Any such kernel whose eigenvalues are all $1$ is called an elementary DPP. According to (Hough et al., 2006, Theorem 7), a DPP with marginal kernel $K$ is a mixture of all $2^{N}$ possible elementary DPPs:

This perspective leads to an efficient DPP sampling algorithm, where a set $J$ is first chosen according to its mixture weight in Equation (7), and then a simple algorithm is used to sample from $P^{V^{J}}$ (Kulesza and Taskar, 2012, Algorithm 1). In this sense, the index set $J$ is an intermediate hidden variable in the process for generating a sample $Y$.

We can exploit this hidden variable $J$ to develop an EM algorithm for learning $K$. Re-writing the data log-likelihood to make the hidden variable explicit:

These equations follow directly from Equations (6) and (7).

3 Lower bounding the objective

We now introduce an auxiliary distribution, $q(J\mid Y_{i})$, and deploy it with Jensen’s inequality to lower-bound the likelihood objective. This is a standard technique for developing EM schemes for dealing with hidden variables Neal and Hinton (1998). Proceeding in this direction:

The function $F(q,V,\Lambda)$ can be expressed in either of the following two forms:

where $H$ is entropy. Consider optimizing this new objective by coordinate ascent. From Equation (11) it is clear that, holding $V,\Lambda$ constant, $F$ is concave in $q$. This follows from the concavity of $\mathbf{KL}$ divergence. Holding $q$ constant in Equation (12) yields the following function:

This expression is concave in $\lambda_{j}$, since $\log$ is concave. However, it is not concave in $V$ due to the non-convex $V^{\top}V=I$ constraint. We describe in Section 3.5 one way to handle this.

To summarize, coordinate ascent on $F(q,V,\Lambda)$ alternates the following “expectation” and “maximization” steps; the first is concave in $q$, and the second is concave in the eigenvalues:

4 E-step

The E-step is easily solved by setting $q(J\mid Y_{i})=p_{K}(J\mid Y_{i})$, which minimizes the KL divergence. Interestingly, we can show that this distribution is itself a conditional DPP, and hence can be compactly described by an $N\times N$ kernel matrix. Thus, to complete the E-step, we simply need to construct this kernel. Lemma 17 (see Appendix A for a proof) gives an explicit formula. Note that $q$’s probability mass is restricted to sets of a particular size $k$, and hence we call it a $k$-DPP. A $k$-DPP is a variant of DPP that can also be efficiently sampled from and marginalized, via modifications of the standard DPP algorithms. (See Appendix A and Kulesza and Taskar (2011b) for more on $k$-DPPs.)

At the completion of the E-step, $q(J\mid Y_{i})$ with $|Y_{i}|=k$ is a $k$-DPP with (non-marginal) kernel $Q^{Y_{i}}$:

5 M-step

The M-step update for the eigenvalues is a closed-form expression with no need for projection. Taking the derivative of Equation (13) with respect to $\lambda_{j}$, setting it equal to zero, and solving for $\lambda_{j}$:

The exponential-sized sum here is impractical, but we can eliminate it. Recall from Lemma 17 that $q(J\mid Y_{i})$ is a $k$-DPP with kernel $Q^{Y_{i}}$. Thus, we can use $k$-DPP marginalization algorithms to efficiently compute the sum over $J$. More concretely, let $\hat{V}$ represent the eigenvectors of $Q^{Y_{i}}$, with $\hat{v}_{r}(j)$ indicating the $j$th element of the $r$th eigenvector. Then the marginals are:

which allows us to compute the eigenvalue updates in time $O(nNk^{2})$, for $k=\max_{i}|Y_{i}|$. (See Appendix B for the derivation of Equation (19) and its computational complexity.) Note that this update is self-normalizing, so explicit enforcement of the $0\leq\lambda_{j}\leq 1$ constraint is unnecessary. There is one small caveat: the $Q^{Y_{i}}$ matrix will be infinite if any $\lambda_{j}$ is exactly equal to 1 (due to $R$ in Equation (17)). In practice, we simply tighten the constraint on $\boldsymbol{\lambda}$ to keep it slightly below 1.

Turning now to the M-step update for the eigenvectors, the derivative of Equation (13) with respect to $V$ involves an exponential-size sum over $J$ similar to that of the eigenvalue derivative. However, the terms of the sum in this case depend on $V$ as well as on $q(J\mid Y_{i})$, making it hard to simplify. Yet, for the particular case of the initial gradient, where we have $q=p$, simplification is possible:

where $H^{Y_{i}}$ is the $|Y_{i}|\times|Y_{i}|$ matrix $V_{Y_{i}}R^{2}V_{Y_{i}}^{\top}$ and $V_{Y_{i}}=(U^{Y_{i}})^{\top}$. $B_{Y_{i}}$ is a $N\times|Y_{i}|$ matrix containing the columns of the $N\times N$ identity corresponding to items in $Y_{i}$; $B_{Y_{i}}$ simply serves to map the gradients with respect to $V_{Y_{i}}$ into the proper positions in $V$. This formula allows us to compute the eigenvector derivatives in time $O(nNk^{2})$, where again $k=\max_{i}|Y_{i}|$. (See Appendix C for the derivation of Equation (20) and its computational complexity.)

Equation (20) is only valid for the first gradient step, so in practice we do not bother to fully optimize $V$ in each M-step; we simply take a single gradient step on $V$. Ideally we would repeatedly evaluate the M-step objective, Equation (13), with various step sizes to find the optimal one. However, the M-step objective is intractable to evaluate exactly, as it is an expectation with respect to an exponential-size distribution. In practice, we solve this issue by performing an E-step for each trial step size. That is, we update $q$’s distribution to match the updated $V$ and $\Lambda$ that define $p_{K}$, and then determine if the current step size is good by checking for improvement in the likelihood $\mathcal{L}$.

There is also the issue of enforcing the non-convex constraint $V^{\top}V=I$. We could project $V$ to ensure this constraint, but, as previously discussed for eigenvalues, projection steps often lead to poor local optima. Thankfully, for the particular constraint associated with $V$, more sophisticated update techniques exist: the constraint $V^{\top}V=I$ corresponds to optimization over a Stiefel manifold, so the algorithm from (Edelman et al., 1998, Page 326) can be employed. In practice, we simplify this algorithm by negelecting second-order information (the Hessian) and using the fact that the $V$ in our application is full-rank. With these simplifications, the following multiplicative update is all that is needed:

where $\exp$ denotes the matrix exponential and $\eta$ is the step size. Algorithm 2 summarizes the overall EM method. As previously mentioned, assuming $T_{1}$ iterations until convergence and an average of $T_{2}$ iterations to find a step size, its overall runtime is $O(T_{1}nNk^{2}+T_{1}T_{2}N^{3})$. The first term in this complexity comes from the eigenvalue updates, Equation (19), and the eigenvector derivative computation, Equation (20). The second term comes from repeatedly computing the Stiefel manifold update of $V$, Equation (21), during the step size search.

Experiments

For the second initialization type, we employ a form of moment matching. Let $m_{i}$ and $m_{ij}$ represent the normalized frequencies of single items and pairs of items in the training data:

Consider a product recommendation task, where the ground set comprises $N$ products that can be added to a particular category (e.g., toys or safety) in a baby registry. A very simple recommendation system might suggest products that are popular with other consumers; however, this does not account for negative interactions: if a consumer has already chosen a carseat, they most likely will not choose an additional carseat, no matter how popular it is with other consumers. DPPs are ideal for capturing such negative interactions. A learned DPP could be used to populate an initial, basic registry, as well as to provide live updates of product recommendations as a consumer builds their registry.

To test our DPP learning algorithms, we collected a dataset consisting of $29{,}632$ baby registries from Amazon.com, filtering out those listing fewer than $5$ or more than $100$ products. Amazon characterizes each product in a baby registry as belonging to one of $18$ categories, such as “toys” and“safety”. For each registry, we created sub-registries by splitting it according to these categories. (A registry with $5$ toy items and $10$ safety items produces two sub-registries.) For each category, we then filtered down to its top $100$ most frequent items, and removed any product that did not occur in at least $100$ sub-registries. We discarded categories with $N<25$ or fewer than $2N$ remaining (non-empty) sub-registries for training. The resulting $13$ categories have an average inventory size of $N=71$ products and an average number of sub-registries $n=8{,}585$. We used $70$% of the data for training and $30$% for testing. Note that categories such as “carseats” contain more diverse items than just their namesake; for instance, “carseats” also contains items such as seat back kick protectors and rear-facing baby view mirrors. See Appendix D for more dataset details and for quartile numbers for all of the experiments.

If instead of the Wishart initialization we use the moments-matching initializer, this alleviates KA’s projection problem, as it provides a starting point closer to the true kernel. With this initializer, KA and EM have comparable test log-likelihoods (average EM gain of 0.4%). However, the moments-matching initializer is not a perfect fix for the KA algorithm in all settings. For instance, consider a data-poor setting, where for each category we have only $n=2N$ training examples. In this case, even with the moments-matching initializer EM has a significant edge over KA, as shown in Figure 1b: EM gains an average of $4.5$%, with a maximum gain of $16.5$% for the safety category.

To give a concrete example of the advantages of EM training, Figure 2a shows a greedy approximation (Nemhauser et al., 1978, Section 4) to the most-likely ten-item registry in the category “safety”, according to a Wishart-initialized EM model. The corresponding KA selection differs from Figure 2a in that it replaces the lens filters and the head support with two additional baby monitors: “Motorola MBP36 Remote Wireless Video Baby Monitor”, and “Summer Infant Baby Touch Digital Color Video Monitor”. It seems unlikely that many consumers would select three different brands of video monitor.

Having established that EM is more robust than KA, we conclude with an analysis of runtimes. Figure 2b shows the ratio of KA’s runtime to EM’s for each category. As discussed earlier, EM is asymptotically faster than KA, and we see this borne out in practice even for the moderate values of $N$ and $n$ that occur in our registries dataset: on average, EM is $2.1$ times faster than KA.

Conclusion

We have explored learning DPPs in a setting where the kernel $K$ is not assumed to have fixed values or a restrictive parametric form. By exploiting $K$’s eigendecomposition, we were able to develop a novel EM learning algorithm. On a product recommendation task, we have shown EM to be faster and more robust than the naive approach of maximizing likelihood by projected gradient. In other applications for which modeling negative interactions between items is important, we anticipate that EM will similarly have a significant advantage.

This work was supported in part by ONR Grant N00014-10-1-0746.

Appendix A Proof of Lemma 1

Lemma 17 gives the exact form of $q$’s kernel. Before giving the proof, we briefly note that $q$ differs slightly from the typical DPPs we have seen thus far, in its conditional nature. More precisely, for a set $Y_{i}$ of size $k$, $q$ qualifies as a $k$-DPP, a DPP conditioned on sets of size $k$. Formally, a $k$-DPP with (non-marginal) kernel $L$ assigns probability $\propto\det(L_{Y})$ for $|Y|=k$, and probability zero for $|Y|\neq k$. As for regular DPPs, a $k$-DPP can be efficiently sampled from and marginalized, via modifications of the standard DPP algorithms. For example, the normalization constant for a $k$-DPP is given by the identity $\sum_{Y:|Y|=k}\det(L_{Y})=e_{k}^{N}(L)$, where $e_{k}^{N}(L)$ represents the $k$th-order elementary symmetric polynomial on the eigenvalues of $L$ Kulesza and Taskar (2011b). Baker and Harwell (1996)’s “summation algorithm” computes $e_{k}^{N}(L)$ in $O(Nk)$ time. In short $k$-DPPs enjoy many of the advantages of DPPs. Their identical parameterization, in terms of a single kernel, makes our E-step simple, and their normalization and marginalization properties are useful for the M-step updates.

Since the E-step is an unconstrained KL divergence minimization, we have:

where the proportionality follows because $Y_{i}$ is held constant in the conditional $q$ distribution. Recalling Equation (9), notice that $p_{K}(Y_{i}\mid J)$ can be re-expressed as follows:

This follows from the identity $\det(AA^{\top})=\det(A^{\top}A)$, for any full-rank square matrix $A$. The subsequent swapping of $J$ and $Y_{i}$, once $V^{\top}$ is re-written as $U$, does not change the indexed submatrix.

where $P^{U^{Y_{i}}}$ represents an elementary DPP, just as in Equation (6), but over $J$ rather than $Y$. Multiplying this expression by a term that is constant for all $J$ maintains proportionality and allows us to simplify the the $p_{K}(J)$ term. Taking the definition of $p_{K}(J)$ from Equation (9):

Having eliminated all dependence on $j\notin J$, it is now possible to express $q(J\mid Y_{i})$ as the $J$ principal minor of a PSD kernel matrix (see $Q^{Y_{i}}$ in the statement of the lemma). Thus, $q$ is a $k$-DPP. ∎

Appendix B M-Step eigenvalue updates

We can exploit standard $k$-DPP marginalization formulas to efficiently compute the eigenvalue updates for EM. Specifically, the exponential-size sum over $J$ from Equation (18) can be reduced to the computation of an eigendecomposition and several elementary symmetric polynomials on the resulting eigenvalues. Let $e_{k-1}^{-j}(Q^{Y_{i}})$ be the $(k-1)$-order elementary symmetric polynomial over all eigenvalues of $Q^{Y_{i}}$ except for the $j$th one. Then, by direct application of (Kulesza and Taskar, 2012, Equation 205), $q$’s singleton marginals are:

As previously noted, elementary symmetric polynomials can be efficiently computed using Baker and Harwell (1996)’s “summation algorithm”.

We can further reduce the complexity of this formula by noting that rank of the $N\times N$ matrix $Q^{Y_{i}}=RZ^{Y_{i}}R$ is at most $|Y_{i}|$. Because $Q^{Y_{i}}$ only has $|Y_{i}|$ non-zero eigenvalues, it is the case that, for all $r$:

Recalling that the eigenvectors and eigenvalues of $Q^{Y_{i}}$ are denoted $\hat{V},\hat{\Lambda}$, the computation of the singleton marginals of $q$ that are necessary for the M-step eigenvalue updates can be written as follows:

Getting $\hat{V}$ via Equation (32) is an $O(N|Y_{i}|^{2})$ operation, given the eigendecomposition of $H^{Y_{i}}$. Since this eigendecomposition is an $O(|Y_{i}|^{3})$ operation, it is dominated by the $O(N|Y_{i}|^{2})$. To compute Equation (31) for all $j$ and requires only $O(Nk)$ time, given $\hat{V}$. Thus, letting $k=\max_{i}|Y_{i}|$, the size of the largest example set, the overall complexity of the eigenvalue updates is $O(nNk^{2})$.

Appendix C M-Step eigenvector gradient

The $p_{K}(J)$ term does not depend on the eigenvectors, so we only have to be concerned with the $p_{K}(Y_{i}\mid J)$ term when computing the eigenvector derivatives. Recall that this term is defined as follows:

Applying standard matrix derivative rules such as (Petersen and Pedersen, 2012, Equation 55), the gradient of the M-step objective with respect to entry $(a,b)$ of $V$ is:

Recall that $Z^{Y_{i}}$ is defined to be $U^{Y_{i}}(U^{Y_{i}})^{\top}$, where $U=V^{\top}$. The $p_{K}(Y_{i}\mid J)$ portion of the M-step objective, Equation (34), can be re-written in terms of $Z^{Y_{i}}$:

Taking the gradient of the M-step objective with respect to $Z^{Y_{i}}$:

Plugging in the $k$-DPP form of $q(J\mid Y_{i})$ derived in the main body of the paper:

Recall from the background section the identity used to normalize a $k$-DPP, and consider taking its derivative with respect to $Z^{Y_{i}}$:

Note that this relationship is only true at the start of the M-step, before $V$ (and hence $Z$) undergoes any gradient updates; a gradient step for $V$ would mean that $Q^{Y_{i}}$, which remains fixed during the M-step, could no longer can be expressed as $RZ^{Y_{i}}R$. Thus, the formula we develop in this section is only valid for the first gradient step.

Plugging Equation (40) back into Equation (39):

Multiplying this by the derivative of $Z^{Y_{i}}$ with respect to $V$ and summing over $i$ gives the final form of the gradient with respect to $V$. Thus, we can compute the value of the first gradient on $V$ exactly in polynomial time.

C.2 Faster computation of the first gradient

Recall from Appendix B that the rank of the $N\times N$ matrix $Q^{Y_{i}}=RZ^{Y_{i}}R$ is at most $|Y_{i}|$ and that its non-zero eigenvalues are identical to those of the $|Y_{i}|\times|Y_{i}|$ matrix $H^{Y_{i}}=V_{Y_{i}}R^{2}V_{Y_{i}}^{\top}$. Since the elementary symmetric polynomial $e_{k}^{N}$ depends only on the eigenvalues of its argument, this means $H^{Y_{i}}$ can substitute for $Q^{Y_{i}}$ in Equation (41), if we change variables back from $Z$ to $V$:

where the $i$-th term in the sum is assumed to index into the $Y_{i}$ rows of the $V$ derivative. Further, because $H$ is only size $|Y_{i}|\times|Y_{i}|$:

Plugging this back into Equation (42) and applying standard matrix derivative rules:

Thus, the initial M-step derivative with respect to $V$ can be more efficiently computed via the above equation. Specifically, the matrix $H^{Y_{i}}$ can be computed in time $O(N|Y_{i}|^{2})$, since $R$ is a diagonal matrix. It can be inverted in time $O(|Y_{i}|^{3})$, which is dominated by $O(N|Y_{i}|^{2})$. Thus, letting $k=\max_{i}|Y_{i}|$, the size of the largest example set, the overall complexity of computing the eigenvector gradient in Equation (44) is $O(nNk^{2})$.

Appendix D Baby registry details

Figure 3a and Figure 3b contain details, referred to in the main body of the paper paper, about the baby registry dataset and the learning methods’ performance on it.

References