Title: EFFICIENT TRAINING ON VERY LARGE CORPORA VIA GRAMIAN ESTIMATION

Abstract: We study the problem of learning similarity functions over very large corpora using neural network embedding models. These models are typically trained using SGD with random sampling of unobserved pairs, with a sample size that grows quadratically with the corpus size, making it expensive to scale. We propose new efficient methods to train such models without sampling unobserved pairs. Inspired by matrix factorization, our approach relies on adding a global quadratic penalty and expressing this term as the inner-product of two Gramians. We show that the gradient of this term can be efficiently computed by maintaining estimates of the Gramians, and develop variance reduction schemes to improve the quality of the estimates. We conduct large-scale experiments that show a significant improvement in training time and generalization performance compared to sampling methods.

1. INTRODUCTION:

We consider the problem of learning a similarity function h : X × Y → R, that maps each pair of items, represented by their feature vectors (x, y) ∈ X × Y, to a real number h(x, y), representing their similarity. We will refer to x and y as the left and right feature vectors, respectively. Many problems can be cast in this form: In natural language processing, x represents a context (e.g. a bag of words), y represents a candidate word, and the target similarity measures the likelihood to observe y in context x (Mikolov et al., 2013; Pennington et al., 2014; Levy & Goldberg, 2014) . In recommender systems, x represents a user query, y represents a candidate item, and the target similarity is a measure of relevance of item y to query x, e.g. a movie rating (Agarwal & Chen, 2009) , or the likelihood to watch a given movie (Hu et al., 2008; Rendle, 2010) . Other applications include image similarity, where x and y are pixel-representations of images (Bromley et al., 1993; Chechik et al., 2010; Schroff et al., 2015) , and network embedding models (Grover & Leskovec, 2016; Qiu et al., 2018) , where x and y are nodes in a graph and the similarity is whether an edge connects them.

A popular approach to learning similarity functions is to train an embedding representation of each item, such that items with high similarity are mapped to vectors that are close in the embedding space. A common property of such problems is that only a small subset of all possible pairs X × Y is present in the training set, and those examples typically have high similarity. Training exclusively on observed examples has been demonstrated to yield poor generalization performance. Intuitively, when trained only on observed pairs, the model places the embedding of a given item close to similar items, but does not learn to place it far from dissimilar ones (Shazeer et al., 2016; Xin et al., 2017) . Taking into account unobserved pairs is known to improve the embedding quality in many applications, including recommendation (Hu et al., 2008; Yu et al., 2017) and word analogy tasks (Shazeer et al., 2016) . This is often achieved by adding a low-similarity prior on all pairs, which acts as a repulsive force between all embeddings. But because it involves a number of terms quadratic in the corpus size, this term is computationally intractable (except in the linear case), and it is typically optimized using sampling: for each observed pair in the training set, a set of random unobserved pairs is sampled and used to compute an estimate of the repulsive term. But as the corpus size increases, the quality of the estimates deteriorates unless the sample size is increased, which limits scalability.

In this paper, we address this issue by developing new methods to efficiently estimate the repulsive term, without sampling unobserved pairs. Our approach is inspired by matrix factorization models, which correspond to the special case of linear embedding functions. They are typically trained using alternating least squares (Hu et al., 2008) , or coordinate descent methods (Bayer et al., 2017) , which circumvent the computational burden of the repulsive term by writing it as a matrix-inner-product of two Gramians, and computing the left Gramian before optimizing over the right embeddings, and viceversa. Unfortunately, in non-linear embedding models, each update of the model parameters induces a simulateneous change in all embeddings, making it impractical to recompute the Gramians at each iteration. As a result, the Gramian formulation has been largely ignored in the non-linear setting, where models are instead trained using stochastic gradient methods with sampling of unobserved pairs, see Chen et al. (2016) . Vincent et al. (2015) were, to our knowledge, the first to attempt leveraging the Gramian formulation in the non-linear case. They consider a model where only one of the embedding functions is non-linear, and show that the gradient can be computed efficiently in that case. Their result is remarkable in that it allows exact gradient computation, but this unfortunately does not generalize to the case where both embedding functions are non-linear.

Contributions:

We propose new methods that leverage the Gramian formulation in the non-linear case, and that, unlike previous approaches, are efficient even when both left and right embeddings are non-linear. Our methods operate by maintaining stochastic estimates of the Gram matrices, and using different variance reduction schemes to improve the quality of the estimates. We perform several experiments that show these methods scale far better than traditional sampling approaches on very large corpora.

We start by reviewing preliminaries in Section 2, then derive the Gramian-based methods and analyze them in Section 3. We conduct large-scale experiments on the Wikipedia dataset in Section 4, and provide additional experiments in the appendix. All the proofs are deferred to Appendix A.

2. PRELIMINARIES:



2.1. NOTATION AND PROBLEM FORMULATION:

We consider models that consist of two embedding functions u : R d ×X → R k and v : R d ×Y → R k , which map a parameter vector1 θ ∈ R d and feature vectors x, y to embeddings u(θ, x), v(θ, y) ∈ R k . The output of the model is the dot product2 of the embeddings h θ (x, y) = u(θ, x), v(θ, y) , where

•, • denotes the usual inner-product on R k . Low-rank matrix factorization is a special case, in which the left and right embedding functions are linear in x and y. Figure 1 illustrates a non-linear model, in which each embedding function is given by a feed-forward neural network. 3We denote the training set by T = {(x i , y i , s i ) ∈ X × Y × R} i∈{1,...,n} , where x i , y i are the feature vectors and s i is the target similarity for example i. To make notation more compact, we will use u i (θ), v i (θ) as a shorthand for u(θ, x i ), v(θ, y i ), respectively. As discussed in the introduction, we also assume that we are given a low-similarity prior p ij ∈ R for all pairs (i, j) ∈ {1, . . . , n} 2 . Given a differentiable scalar loss function : R × R → R, the objective function is given by

EQUATION

where the first term measures the loss on observed data, the second term penalizes deviations from the prior, and λ is a positive hyper-parameter that trades-off the two terms. To simplify the discussion, we will assume a uniform zero prior p ij as in (Hu et al., 2008) , the general case is treated in Appendix B.

To optimize this objective, existing methods rely on sampling to approximate the second term, and are usually referred to as negative sampling or candidate sampling, see Chen et al. (2016) ; Yu et al. (2017) for a survey. Due to the double sum in (1), the quality of the sampling estimates degrades as the corpus size increases, which can significantly increase training times. This can be alleviated by increasing the sample size, but does not scale to very large corpora.

x

∈ X u(θ, x) ∈ R k . . . . . . . . . . . . u(θ, •) y ∈ Y v(θ, x) ∈ R k . . . . . . . . . . . . v(θ, •) • h θ (x, y) = u(θ, x), v(θ, y)

Figure 1 : A dot-product embedding model for a similarity function on X × Y.

2.2. GRAMIAN FORMULATION:

A different approach to solving (1), widely popular in matrix factorization, is to rewrite the double sum as the inner product of two Gram matrices. Let us denote by U θ ∈ R n×k the matrix of all left embeddings such that u i (θ) is the i-th row of U θ , and similarly for V θ ∈ R n×k . Then denoting the matrix inner-product by A, B = i,j A ij B ij , we can rewrite the double sum in (1) as:

EQUATION

Now, using the adjoint property of the inner product, we have

U θ V θ , U θ V θ = U θ U θ , V θ V θ ,

and if we denote by u ⊗ u the outer product of a vector u by itself, and define the Gram matrices

EQUATION

the penalty term becomes

g(θ) = G u (θ), G v (θ) .

(4) The Gramians are k × k PSD matrices, where k, the dimension of the embedding space, is much smaller than n -typically k is smaller than 1000, while n can be arbitrarily large. Thus, the Gramian formulation (4) has a much lower computational complexity than the double sum formulation (2), and this reformulation is at the core of alternating least squares and coordinate descent methods (Hu et al., 2008; Bayer et al., 2017) , which operate by computing the exact Gramian for one side, and solving for the embeddings on the other. However, these methods do not apply in the non-linear setting due to the implicit dependence on θ, as a change in the model parameters simultaneously changes all embeddings on both sides, making it intractable to recompute the Gramians at each iteration, so the Gramian formulation has not been used when training non-linear models. In the next section, we show that it can in fact be leveraged in the non-linear case.

3. TRAINING EMBEDDING MODELS USING GRAMIAN ESTIMATES:

In order to leverage the Gramian formulation in the non-linear case, we start by rewriting the objective function (1) in terms of the Gramians defined in (3). Let

EQUATION

then (1) is equivalent to

min θ∈R d 1 n n i=1 [f i (θ) + λg i (θ)]

. Intuitively, for each example i, -∇f i (θ) pulls the embeddings u i and v i close to each other (assuming a high similarity s i ), while -∇g i (θ) creates a repulsive force between u i and all embeddings {v j } j∈{1,...,n} , and between v i and all {u j } j∈{1,...,n} , see Appendix D for further discussion, and illustration of the effect of this term. While the Gramians are expensive to recompute at each iteration, we can maintain PSD estimates Ĝu , Ĝv of G u (θ), G v (θ). Then the gradient of g(θ) can be approximated by the gradient (w.r.t. θ) of ĝi (θ, Ĝu , Ĝv ) :

EQUATION

as stated in the following proposition.

Proposition 1. If i is drawn uniformly in {1, . . . , n}, and Ĝu , Ĝv are unbiased estimates of G u (θ), G v (θ) and independent of i, then ∇ θ ĝi (θ, Ĝu , Ĝv ) is an unbiased estimate of ∇g(θ).

In a mini-batch setting, one can further average ĝi over a batch of examples i ∈ B (which we do in our experiments), but we will omit batches to keep the notation concise. Next, we propose several methods for computing Gramian estimates Ĝu , Ĝv , and discuss their tradeoffs. Since each Gramian can be written as a sum of rank-one terms, e.g. G u (θ) = 1 n n i=1 u i (θ) ⊗ u i (θ), a simple unbiased estimate can be obtained by sampling one term (or a batch) from this sum. We improve on this by using different variance reduction methods, which we discuss in the next two sections.

3.1. STOCHASTIC AVERAGE GRAMIAN:

Our first method is inspired by the stochastic average gradient (SAG) method (Roux et al., 2012; Defazio et al., 2014; Schmidt et al., 2017) , which reduces the variance of the gradient estimates by maintaining a cache of individual gradients, and estimating the full gradient using this cache. Since each Gramian is a sum of outer-products (see equation ( 3)), we can apply the same technique to estimate Gramians. For all i ∈ {1, . . . , n}, let ûi , vi be a cache of the left and right embeddings respectively. We will denote by a superscript (t) the value of a variable at iteration t. Let

Ŝ(t) u = 1 n n i=1 û(t) i ⊗ û(t)

i , which corresponds to the Gramian computed with the current caches. At each iteration t, an example i is drawn uniformly at random and the estimate of the Gramian is given by

EQUATION

and similarly for Ĝ(t) v . This is summarized in Algorithm 1, where the model parameters are updated using SGD (line 10), but this update can be replaced with any first-order method. Here β can take one of the following values: β = 1 n , following SAG (Schmidt et al., 2017) , or β = 1, following SAGA (Defazio et al., 2014) . The choice of β comes with trade-offs that we briefly discuss below. We will denote the cone of positive semi-definite k × k matrices by

S k + . Proposition 2. Suppose β = 1 n in (8). Then for all t, Ĝ(t) u , Ĝ(t) v ∈ S k + . Proposition 3. Suppose β = 1 in (8). Then for all t, Ĝ(t) u is an unbiased estimate of G u (θ (t) ).

While taking β = 1 gives an unbiased estimate, note that it does not guarantee that the estimates remain in S k + . In practice, this can cause numerical issues, but can be avoided by projecting Ĝu , Ĝv on S k + , using their eigenvalue decompositions. The per-iteration cost of maintaining the Gramian estimates is O(k) to update the caches, O(k 2 ) to update the estimates Ŝu , Ŝv , Ĝu , Ĝv , and O(k 3 ) for projecting on S k + . Given the small size of the embedding dimension k, O(k 3 ) remains tractable. The memory cost is O(nk), since each embedding needs to be cached (plus a negligible O(k 2 ) for storing the Gramian estimates). This makes SAGram much less expensive than applying the original SAG(A) methods, which require maintaining caches of the gradients, this would incur a O(nd) memory cost, where d is the number of parameters of the model, and can be orders of magnitude larger than the embedding dimension k. However, O(nk) can still be prohibitively expensive when n is very large. In the next section, we propose a different method which does not incur this additional memory cost.

Algorithm 1 SAGram (Stochastic Average Gramian) 1: Input: Training data {(x i , y i , s i )} i∈{1,...,n} , learning rate η > 0. 2: Initialization phase 3: draw θ randomly 4: ûi ← u i (θ), vi ← v i (θ) ∀i ∈ {1, . . . , n} 5: Ŝu ← 1 n n i=1 ûi ⊗ ûi , Ŝv ← 1 n n i=1 vi ⊗ vi 6: repeat 7: Update Gramian estimates (i ∼ Uniform(n)) 8: Ĝu ← Ŝu + β[u i (θ) ⊗ u i (θ) -ûi ⊗ ûi ], Ĝv ← Ŝv + β[v i (θ) ⊗ v i (θ) -vi ⊗ vi ] 9: Update model parameters then update caches (i ∼ Uniform(n)) 10: θ ← θ -η∇ θ [f i (θ) + λĝ i (θ, Ĝu , Ĝv )] 11: Ŝu ← Ŝu + 1 n [u i (θ) ⊗ u i (θ) -ûi ⊗ ûi ], Ŝv ← Ŝv + 1 n [v i (θ) ⊗ v i (θ) -vi ⊗ vi ] 12: ûi ← u i (θ), vi ← v i (θ) 13: until stopping criterion

3.2. STOCHASTIC ONLINE GRAMIAN:

To derive the second method, we reformulate problem (1) as a two-player game. The first player optimizes over the parameters of the model θ, the second player optimizes over the Gramian estimates Ĝu , Ĝv ∈ S k + , and they seek to minimize the respective losses

EQUATION

where ĝi is defined in ( 7), and • F denotes the Frobenius norm. To justify this reformulation, we can characterize its first-order stationary points, as follows.

Proposition 4. (θ, Ĝu , Ĝv ) ∈ R d × S k + × S k

+ is a first-order stationary point for ( 9) if and only if θ is a first-order stationary point for problem (1) and Ĝu

= G u (θ), Ĝv = G v (θ).

Several stochastic first-order dynamics can be applied to problem (9), and Algorithm 2 gives a simple instance where each player implements SGD with a constant learning rate. In this case, the updates of the Gramian estimates (line 7) have a particularly simple form, since

EQUATION

and similarly for Ĝv . One advantage of this form is that each update performs a convex combination between the current estimate and a rank-1 PSD matrix, thus guaranteeing that the estimates remain in S k + , without the need to project. The per-iteration cost of updating the estimates is O(k 2 ), and the memory cost is O(k 2 ) for storing the Gramians, which are both negligible.

Algorithm 2 SOGram (Stochastic Online Gramian) 1: Input: Training data {(x i , y i , s i )} i∈{1,...,n} , learning rates η > 0, α ∈ (0, 1). Update Gramian estimates (i ∼ Uniform)

6: Ĝu ← (1 -α) Ĝu + αu i (θ) ⊗ u i (θ), Ĝv ← (1 -α) Ĝv + αv i (θ) ⊗ v i (θ) 7: Update model parameters (i ∼ Uniform) 8: θ ← θ -η∇ θ [f i (θ) + λĝ i (θ, Ĝu , Ĝv )] 9: until stopping criterion

The update (10) can also be interpreted as computing an online estimate of the Gramian by averaging rank-1 terms with decaying weights, thus we call the method Stochastic Online Gramian. Indeed, we have by induction on t, Ĝ(t)

u = t τ =1 α(1 -α) t-τ u iτ (θ (τ ) ) ⊗ u iτ (θ (τ )

). Intuitively, averaging reduces the variance of the estimator but introduces a bias, and the choice of the hyper-parameter α trades-off bias and variance. The next proposition quantifies this tradeoff under mild assumptions.

Proposition 5. Let Ḡ(t) u = t τ =1 α(1 -α) t-τ G u (θ (τ )

). Suppose that there exist σ, δ > 0 such that for all t,

EQUATION

)

The first assumption simply bounds the variance of single-point estimates, while the second bounds the distance between two consecutive Gramians, a reasonable assumption, since in practice the changes in Gramians vanish as the trajectory θ (τ ) converges. In the limiting case α = 1, Ĝ(t) u reduces to a single-point estimate, in which case the bias (12) vanishes and the variance ( 11) is maximal, while smaller values of α decrease variance and increase bias. This is confirmed in our experiments, as discussed in Section 4.2.

3.3. COMPARISON WITH EXISTING STOCHASTIC METHODS:

We conclude this section by showing that candidate sampling methods (see Chen et al. (2016) ; Yu et al. (2017) for recent surveys) can be reinterpreted in terms of the Gramian formulation (4). These methods work by approximating the double-sum in (1) using a random sample of pairs. Suppose a batch of pairs (i, j) ∈ B × B is sampled4 , and the double sum is approximated by

EQUATION

)

where µ i , ν j are the inverse probabilities of sampling i, j respectively (to guarantee that the estimate is unbiased). Then applying a similar transformation to Section 2.2, one can show that

EQUATION

)

which is equivalent to computing two batch-estimates of the Gramians. Implementing existing methods using ( 14) rather than (13) can decrease their computional complexity in the large batch regime, for the following reason: the double-sum formulation (13) involves a sum of |B||B | dot products of vectors in R k , thus computing its gradient costs O(k|B||B |). On the other hand, the Gramian formulation ( 14) is the inner product of two k × k matrices, each involving a sum over the batch, thus computing its gradient costs O(k 2 max(|B|, |B |)), which is cheaper when the batch size is larger than the embedding dimension k, a common situation in practice. With this formulation, the advantage of SOGram and SAGram becomes clear, as they use more embeddings to estimate Gramians (by caching or online averaging) than would be possible using candidate sampling.

4. EXPERIMENTS:

In this section, we conduct large-scale experiments on the Wikipedia dataset (Wikimedia Foundation). Additional experiments on MovieLens (Harper & Konstan, 2015) are given in Appendix F.

4.1. EXPERIMENTAL SETUP:



Datasets:

We consider the problem of learning the intra-site links between Wikipedia pages. Given a pair of pages (x, y) ∈ X × X , the target similarity is 1 if there is a link from x to y, and 0 otherwise. Here a page is represented by a feature vector x = (x id , x ngrams , x cats ), where x id is (a one-hot encoding of) the page URL, x ngrams is a bag-of-words representation of the set of n-grams of the page's title, and x cats is a bag-of-words representation of the categories the page belongs to. Note that the left and right feature spaces coincide in this case, but the target similarity is not necessarily symmetric (the links are directed edges). We carry out experiments on subsets of the Wikipedia graph corresponding to three languages: Simple English, French, and English, denoted respectively by simple, fr, and en. These subgraphs vary in size, and Table 1 shows some basic statistics for each set. Each set is partitioned into training and validation using a (90%, 10%) split. Models We train non-linear embedding models consisting of a two-tower neural network as in Figure 1 , where the left and right embedding functions map, respectively, the source and destination page features. The two embedding networks have the same structure: the input feature embeddings are concatenated then mapped through two hidden layers with ReLU activations. The input embeddings are shared between the two networks, and their dimensions are 50 for simple, 100 for fr, and 120 for en. The sizes of the hidden layers are [256, 64] for simple and [512, 128] for fr and en.

Training methods The model is trained using a squared error loss, (s, s ) = 1 2 (ss ) 2 , optimized using SAGram, SOGram, and as baseline, SGD with candidate sampling, using different sampling strategies. The experiments reported in this section use a learning rate η = 0.01, a penalty coefficient λ = 10, and batch size 1024. These parameters correspond to the best performance of the baseline methods; we report additional results with different hyper-parameter settings in Appendix E. For SAGram and SOGram, a batch B is used in the Gramian updates (line 8 in Algorithm 1 and line 6 in Algorithm 2, where we use a sum of rank-1 terms over the batch), and another batch B is used in the model parameter update. For the sampling baselines, the double sum is approximated by all pairs in the cross product (i, j) ∈ B × B , and for efficiency, we implement them using the Gramian formulation as discussed in Section 3.3, since we operate in a regime where the batch size is an order of magnitude larger than the embedding dimension k. In the first baseline method, uniform, items are sampled uniformly from the vocabulary (all pages are sampled with the same probability). The other baseline methods implement importance sampling similarly to Bengio & Senecal (2003) ; Mikolov et al. (2013) : in linear, the probability is proportional to the number of occurrences of the page in the training set, and in sqrt, the probability is proportional to the square root of the number of occurrences. Figure 2 : Gramian estimation error on a common trajectory (θ (t) ).

4.2. QUALITY OF GRAMIAN ESTIMATES:

In the first set of experiments, we evaluate the quality of the Gramian estimates using each method.

In order to have a meaningful comparison, we fix a trajectory of model parameters (θ (t) ) t∈{1,...,T } , and evaluate how well each method tracks the true Gramians G u (θ (t) ), G v (θ (t) ) on that common trajectory. This experiment is done on Wikipedia simple (the smallest of the datasets) so that we can compute the exact Gramians by periodically computing the embeddings u i (θ (t) ), v i (θ (t) ) on the full training set at a given time t. We report in Figure 2 the estimation error for each method, measured by the normalized Frobenius distance

Ĝ(t) u -Gu(θ (t) ) F Gu(θ (t) ) F

. In Figure 2a , we can observe that both variants of SAGram yield the best estimates, and that SOGram yields better estimates than the baselines. Among the baseline methods, importance sampling (both linear and sqrt) perform better than uniform. We also vary the batch size to evaluate its impact: increasing the batch size from 128 to 1024 improves the quality of all estimates, as expected, but it is worth noting that the estimates of SOGram with |B| = 128 have comparable quality to baseline estimates with |B| = 1024.

In Appendix E, we show that a similar effect can be observed for gradient estimates, and we make a formal connection between Gramian and gradient estimation errors.

In Figure 2b , we evaluate the bias-variance tradeoff discussed in Section 3.2, by comparing the estimates of SOGram with different learning rates α. We observe that higher values of α suffer from higher variance which persists throughout the trajectory. A lower α reduces the variance but introduces a bias, which is mostly visible during the early iterations.

4.3. IMPACT ON TRAINING SPEED AND GENERALIZATION QUALITY:

In order to evaluate the impact of the Gramian estimation quality on training speed and generalization, we compare the validation performance of SOGram to the sampling baselines, on each dataset (we do not use SAGram due to its prohibitive memory cost for corpus sizes of 1M or more). The models are trained with a fixed time budget of 20 hours for simple, 30 hours for fr and 50 hours for en. We estimate the mean average precision (MAP) at 10, by scoring, every 5 minutes, left items in the validation set against 50K random candidates (exhaustively scoring all candidates is prohibitively expensive, but this gives a reasonable approximation). The results are reported in Figure 3 . Compared to the sampling baselines, SOGram exhibits faster training and better validation performance across all sampling strategies. Table 2 The improvement on simple is modest (between 4% and 10%), which can be explained by the relatively small corpus size (85K unique pages), in which case candidate sampling with a large batch size already yields decent estimates. On the larger corpora, we obtain more significant improvements: between 9% and 15% on fr and between 9% and 19% on en. It's interesting to observe that the best performance is consistently achieved by SOGram with linear importance sampling, even though linear performs slightly worse than other strategies in the baseline. SOGram also has a significant impact on training speed: if we measure the time it takes for SOGram to exceed the final validation performance of each baseline method, this time is a small fraction of the total budget. In our experiments, this fraction is between 10% and 17% for simple, between 23% and 30% for fr, and between 16% and 24% for en. Additional numerical results are provided in Appendix E, where we evaluate the impact of other parameters, such as the effect of the batch size |B|, the learning rate η, and the Gramian learning rate α. For example, we show that the relative improvement of SOGram compared to the baselines is even larger when using smaller batches, and its generalization performance is more robust to the choice of batch size and learning rate.

5. CONCLUSION:

We showed that the Gramian formulation commonly used in low-rank matrix factorization can be leveraged for training non-linear embedding models, by maintaining estimates of the Gram matrices and using them to estimate the gradient. By applying variance reduction techniques to the Gramians, one can improve the quality of the gradient estimates, without relying on large sample size as is done in traditional sampling methods. This leads to a significant impact on training time and generalization quality, as indicated by our experiments. While we focused on problems with very large vocabulary size, where traditional approaches are inefficient, it will be interesting to evaluate our methods on other applications such as word-analogy tasks Mikolov et al. (2013) ; Schnabel et al. (2015) . Another direction of future work is to extend this formulation to a larger family of penalty functions, such as the spherical loss family studied in (Vincent et al., 2015; de Brébisson & Vincent, 2016) .

A PROOFS:

Proof of Proposition 1. Starting from the expression (6) of g(θ) = G u (θ), G v (θ) , and applying the chain rule, we have

EQUATION

where J u (θ) denotes the Jacobian of G u (θ), an order-three tensor given by J u (θ) l,i,j = ∂G u (θ) i,j ∂θ l , l ∈ {1, . . . , d}, i, j ∈ {1, . . . , n},

and J u (θ)[G v (θ)] denotes the vector [ i,j J u (θ) l,i,j G v (θ) i,j ] l∈{1,...,d} . Observing that ĝi (θ, Ĝu , Ĝv ) = Ĝu , u i (θ) ⊗ u i (θ) + Ĝv , v i (θ) ⊗ v i (θ)

, and applying the chain rule, we have

EQUATION

) where J u,i (θ) is the Jacobian of u i (θ) ⊗ u i (θ), and

E i∼Uniform [J u,i (θ)] = 1 n n i=1 J u,i (θ) = J u (θ),

an similarly for J v,i . We conclude by taking expectations in ( 16) and using assumption that Ĝu , Ĝv are independent of i.

Proof of Proposition 2. From ( 8) and the definition of

Ŝ(t) u , we have, Ĝ(t) u = 1 n j =i û(t) j ⊗ û(t) j + 1 n u i (θ (t) ) ⊗ u i (θ (t) ),

which is a sum of matrices in the PSD cone S k + .

Proof of Proposition 3. Denoting by (F t ) t≥0 the filtration generated by the sequence (θ (t) ) t≥0 , and taking conditional expectations in (8), we have

E[ Ĝ(t) u |F t ] = Ŝ(t) u + E i∼Uniform [u i (θ (t) ) ⊗ u i (θ (t) ) - û(t) i ⊗ û(t) i |F t ] = Ŝ(t) u + 1 n n i=1 [u i (θ (t) ) ⊗ u i (θ (t) ) -ûi ⊗ ûi ] = 1 n n i=1 u i (θ (t) ) ⊗ u i (θ (t) ) = G u (θ (t) ). Proof of Proposition 4. (θ, Ĝu , Ĝv ) ∈ R d × S k + × S k

+ is a first-order stationary point of the game if and only if

EQUATION

)

The second and third conditions simply states that ∇ Ĝu L θ 2 ( Ĝu , Ĝv ) and ∇ Ĝv L θ 2 ( Ĝu , Ĝv ) define supporting hyperplanes of S k + at Ĝu , Ĝv , respectively.

Since G u (θ) ∈ S k + , condition ( 18) is equivalent to Ĝu = G u (θ) (and similarly, ( 19) is equivalent to Ĝv = G v (θ)). Using the expression (15) of ∇g, we get that (17-19) is equivalent to ∇f (θ) + λ∇g(θ) = 0.

Published as a conference paper at ICLR 2019 Algorithm 3 SAGram (Stochastic Average Gramian) with low-rank prior 1: Input: Training data {(x i , y i , s i )} i∈{1,...,n} , low-rank priors {q i , r i } i∈{1,...,n} 2: Initialization phase

3: draw θ randomly 4: ûi ← u i (θ), vi ← v i (θ) ∀i ∈ {1, . . . , n} 5: Ŝu ← 1 n n i=1 ûi ⊗ ûi , Ŝv ← 1 n n i=1 vi ⊗ vi 6: Tu ← 1 n n i=1 ûi ⊗ q i , Tv ← 1 n n i=1 vi ⊗ r i 7: repeat 8: Update Gramian estimates (i ∼ Uniform(n)) 9: Ĝu ← Ŝu + β[u i (θ) ⊗ u i (θ) -ûi ⊗ ûi ], Ĝv ← Ŝv + β[v i (θ) ⊗ v i (θ) -vi ⊗ vi ] 10:

Update weighted embedding estimates 11:

Ĥu ← Tu + λ[(u i (θ) -ûi ) ⊗ q i ] 12: Ĥv ← Tv + λ[(v i (θ) -vi ) ⊗ r i ] 13: Update model parameters then update caches (i ∼ Uniform(n)) 14: θ ← θ -η∇ θ [f i (θ) + λpĝ P (θ, Ĝu , Ĝv , Ĥu , Ĥv )] 15: Ŝu ← Ŝu + 1 n [u i (θ) ⊗ u i (θ) -ûi ⊗ ûi ], Ŝv ← Ŝv + 1 n [v i (θ) ⊗ v i (θ) -vi ⊗ vi ]

16:

Tu ← Tu + 1 n [(u i (θ) -ûi ) ⊗ q i ], Tv ← Tv + 1 n [(v i (θ) -vi ) ⊗ r i ] 17: ûi ← u i (θ), vi ← v i (θ) 18: until stopping criterion Algorithm 4 SOGram (Stochastic Online Gramian) with low-rank prior 1: Input: Training data {(x i , y i , s i )} i∈{1,...,n} , low-rank priors {q i , r i } i∈{1,...,n} 2: Initialization phase 3: draw θ randomly 4: Ĝu , Ĝv ← 0 k×k 5: repeat 6: Update Gramian estimates (i ∼ Uniform(n)) 7: Ĝu ← (1 -α) Ĝu + αu i (θ) ⊗ u i (θ), Ĝv ← (1 -α) Ĝv + αv i (θ) ⊗ v i (θ) 8:

Update weighted embedding estimates

9: Ĥu ← (1 -α) Ĥu + αu i (θ) ⊗ q i , Ĥv ← (1 -α) Ĥv + αv i (θ) ⊗ r i 10: Update model parameters (i ∼ Uniform(n)) 11: θ ← θ -η∇ θ [f i (θ) + λpĝ P (θ, Ĝu , Ĝv , Ĥu , Ĥv )] 12: until stopping criterion C NON-UNIFORM WEIGHTS

In additional to using a non-uniform prior, it can also be desirable to use non-uniform weights in the penalty term, for example to balance the contribution of frequent and infrequent items to the penalty term. We discuss how to adapt our algorithms to the non-uniform weights case. Suppose that the penalty function is given by

g W (θ) = 1 n 2 n i=1 n j=1 a i b j u i (θ), v j (θ) 2 ,

where a i , b j are positive left and right weights, respectively. Here we used a superscript W in g W to disambiguate the uniform-weight case. Then using a similar transformation to Section 2.2, we can rewrite g W as follows:

g W (θ) = 1 n n i=1 a i u i ⊗ u i , 1 n n j=1 b j v j ⊗ v j ,

i.e. g W is the inner-product of two weighted Gramians. Both SAGram and SOGram can be generalized to this case, by maintaining estimates of the weighted Gramians, one simply needs to scale the contribution of each term u i ⊗ u i by the appropriate embedding weight a i (and similarly for the right embedding).

Remark Here we discussed the case of a rank-one weight matrix, i.e. when the unobserved weight matrix can be written as W = a ⊗ b for a given left and right weight vectors a, b. The weight matrix cannot be arbitrary (as specifying n 2 individual weights is prohibitively expensive in many applications such as the experiments of this paper), thus one needs a consice description of the weights matrix. One such description is the sum of a sparse and low-rank matrix, and one can generalize SAGram and SOGram to this case: the sparse part of the weight matrix can be optimized explicitly, and the low-rank part can be optimized using weighted Gramians, by generalizing the argument of the previous paragraph.

D INTERPRETATIONS OF THE GRAMIAN PENALTY TERM:

In this section, we briefly discuss different interpretations of the Gramian inner-product g(θ). Starting from the expression (4) of g(θ) and the definition (3) of the Gram matrices, we have

EQUATION

which is a quadratic form in the left embeddings u i (and similarly for v j , by symmetry). In particular, the partial derivative of the Gramian term with respect to an embedding u i is

∂g(θ) ∂u i = 2 n G v (θ)u i (θ) = 2 n   1 n n j=1 v j (θ) ⊗ v j (θ)   u i (θ).

Each term (v j ⊗ v j )u i = v j v j , u i is simply the projection of u i on v j (scaled by v j 2 ). Thus the gradient of g(θ) with respect to u i is an average of scaled projections of u i on each of the right embeddings v j , and moving in the direction of the negative gradient simply moves u i away from regions of the embedding space with a high density of right embeddings. This corresponds to the intuition discussed in the introduction: the purpose of the g(θ) term is precisely to push left and right embeddings away from each other, to avoid placing embeddings of dissimilar items near each other, a phenomenon referred to as folding of the embedding space (Xin et al., 2017) . In order to illustrate the effect of this term on the embedding distributions, we visualize, in Figure 4 , the distribution of the inner product u i (θ (t) ), v j (θ (t) ) , for random pairs (i, j), and for observed pairs (i = j), and how these distributions change as t increases. The plots are generated for the Wikipedia en model described in Section 4, trained with SOGram (α = 0.01), with two different values of the penalty coefficient, λ = 10 -2 and λ = 10. In both cases, the distribution for observed pairs remains concentrated around values close to 1, as one expects (recall that the target similarity is 1 for observed pairs, i.e. pairs of connected pages in the Wikipedia graph). The distributions for random pairs, however, are very different: with λ = 10, the distribution quickly concentrates around a value close to 0, while with λ = 10 -2 the distribution is more flat, and a large proportion of pairs have a high inner-product. This indicates that with a lower λ, the model is more likely to fold, i.e. place embeddings of unrelated items near each other. This is consistent with the validation MAP, reported in Figure 5 . With λ = 10 -2 , the validation MAP increases very slowly, and remains two orders of magnitude smaller than the model trained with λ = 10. The figure also shows that when λ is too large (λ = 10 3 ), the model is over-regularized and the MAP decreases.

To conclude this section, we note that our methods also apply to a related regularizer introduced in (Zhang et al., 2017) , called Global Orthogonal Regularization. The authors argue that when learning feature embedding representations, spreading out the embeddings is helpful for generalization, and propose to match the second moment of each embedding distribution with that of the uniform distribution. Formally, and using our notation, they use the penalty term max(g u (θ), 1/k)+max(g v (θ), 1/k),

where k is the embedding dimension, g u (θ) = 1 n 2 n i=1 n j=1 u i , u j 2 , and similarly for g v . They optimize this term using candidate sampling. We can also apply the same Gramian transformation as in Section 2.2 to write g

u (θ) = G u (θ), G u (θ) , and g v (θ) = G v (θ), G v (θ)

, and we can similarly apply SAGram and SOGram to estimate both Gramians. Formally, the difference here is that one would penalize the inner-product of each Gramian with itself, instead of the inner-product of the two. One advantage of this regularizer is that it applies to a broader class of models, as it does not require the output of the model to be the dot-product of two embedding functions.

E FURTHER EXPERIMENTS ON WIKIPEDIA E.1 QUALITY OF GRADIENT ESTIMATES:

The experiments in Section 4 indicate that our methods give better estimates of the Gramians, and a natural question is how this affects gradient estimation quality. First, one can make a formal connection between the two. Since ∂g(θ) We also evaluate the quality of the Gramian estimates, but due to the large vocabulary size in en, computing the exact Gramians is no longer feasible, so we approximate it using a large sample of 1M embeddings. The results are reported in Figure 8 , which shows the normalized Frobenius distance between the Gramian estimates Ĝu and (the large sample approximation of) the true Gramian G u .

∂u i = 2 n G V (θ)u i (θ) ∂ĝ(θ) ∂u i = 2 n ĜV u i (θ),

The results are similar to the experiment on simple: with a lower α, the estimation error is initially high, but decays to a lower value as training progresses, which can be explained by the bias-variance tradeoff discussed in Proposition 5. The tradeoff is affected by the trajectory of the true Gramians: smaller changes in the Gramians (captured by the parameter δ in Proposition 5) induce a smaller bias. In particular, changing the learning rate η of the main algorithm can affect the performance of the Gramian estimates by affecting the rate of change of the true Gramians. To investiage this effect, we ran the same experiment with two different learning rates, η = 0.01 as in Section 4, and a lower learning rate η = 0.002. The errors converge to similar values in both cases, but the error decay occurs much faster with smaller η, which is consistent with our analysis.

E.4 EFFECT OF BATCH SIZE AND LEARNING RATE:

In this section, we explore the effect of the batch size |B| and learning rate η on the performance of SOGram compared to the baselines. We ran the Wikipedia en experiment with different values of these hyperparameters, and report the final validation MAP in Tables 3 and 4 , which correspond to batch size 128 and 512 respectively. We can make several observations. First, the best performance is consistently achieved by SOGram with learning rate α = 0.001. Second, the relative improvement compared to the baseline is, in general, larger for smaller batch sizes. This can be explained intuitively by the fact that because of online averaging, the quality of the Gramian estimates with SOGram suffers less than with the sampling baseline. Finally, we can also observe that the final performance also seems more robust to the choice of batch size and learning rate, compared to the baseline. For example, with the larger learning rate η = 0.02, the performance degrades for all methods, but the drop in performance for the baseline is much more significant than for the SOGram methods.

F EXPERIMENT ON MOVIELENS DATA:

In this section, we report experiments on a regression task on MovieLens.

Dataset The MovieLens dataset consists of movie ratings given by a set of users. In our notation, the left features x represent a user, the right features y represent an item, and the target similarity is the rating of movie y by user x. The data is partitioned into a training and a validation set using a (80%-20%) split. Model We train a two-tower neural network model, as described in Figure 1 , where each tower consists of an input layer, a hidden layer, and output embedding dimension k = 35. The left tower takes as input a one-hot encoding of a unique user id, and the right tower takes as input one-hot encodings of a unique movie id, the release year of the movie, and a bag-of-words representation of the genres of the movie. These input embeddings are concatenated and used as input to the right tower.

Methods The model is trained using a squared loss (s, s ) = 1 2 (ss ) 2 , using SOGram with different values of α, and sampling as a baseline. We use a learning rate η = 0.05, and penalty coefficient λ = 1. We measure mean average precision on the trainig set and validation set, following the same procedure described in Section 4. The results are given in Figure 9 .

Results:

The results are similar to those reported on the Wikipedia simple dataset, which is comparable in corpus size and number of observations to MovieLens. The best validation mean average precision is achieved by SOGram with α = 0.1 (for an improvement of 2.9% compared to the sampling baseline), despite its poor performance on the training set, which indicates that better 

:

In many applications, it is desirable for the two embedding functions u, v to share certain parameters, e.g. embeddings of categorical features common to left and right items; hence, we use the same θ for both.

This also includes cosine similarity models when the embedding functions u, v are normalized.

We focus on models with this dot-product structure since they allow efficient retrieval: given x ∈ X , finding items y ∈ Y with high similarity to x reduces to a maximum inner-product search problem (MIPS), which can be approximated efficiently(Shrivastava & Li, 2014;Neyshabur & Srebro, 2015).

The simplest variant uses a uniform distribution, but other methods have been proposed, such adaptive sampling(Bengio & Senecal, 2008;Bai et al., 2017), and importance sampling(Bengio & Senecal, 2003).

Figures and Tables:

θ randomly, Ĝu , Ĝv ← 0 k×k 4: repeat 5:

SOGram α = 0. 001 SOGram α = 0. 005 SOGram α = 0. 01 SOGram α = 0. 1 (b) SOGram with different averaging rates.

(a) λ = 10 -2 , observed pairs. (b) λ = 10 -2 , random pairs. (c) λ = 10, observed pairs. (d) λ = 10, random pairs.

Figure 4: Evolution of the inner product distribution u i (θ (t) ), v j (θ (t) ) in the Wikipedia en model trained with different penalty coefficients λ, for observed pairs (left) and random pairs (right).

Figure 5: Mean Average Precision of the Wikipedia en model, trained with different values of the penalty coefficient λ.

Figure 7: Validation MAP of SOGram and linear importance sampling, on Wikipedia en, for different values of the Gramian learning rates α.

Figure 8: Gramian estimation error on en, for SOGram with different values of α, and different learning rates. The left and right figures correspond respectively to η = 0.01 and η = 0.002.

Figure 9: Mean average precision at 10 on the training set (left) and the validation set (right), for different methods, on the MovieLens dataset.

Corpus sizes for each training set.

Table content: <table><tr><td colspan="4">language # pages # links # ngrams</td><td># cats</td></tr><tr><td>simple</td><td>85K</td><td>4.6M</td><td>8.3K</td><td>6.1K</td></tr><tr><td>fr</td><td>1.8M</td><td>142M</td><td colspan="2">167.4K 125.3K</td></tr><tr><td>en</td><td>5.3M</td><td>490M</td><td colspan="2">501.0K 403.4K</td></tr></table>

summarizes the relative improvement of the final validation MAP.

Table content: <table><tr><td/><td>0.032</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>0.18</td><td/><td/><td/></tr><tr><td/><td>0.030</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td/><td/><td>0.11</td><td/><td/><td/><td/><td/><td/><td>0.16</td><td/><td/><td/></tr><tr><td/><td>0.028</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td/><td/><td>0.10</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Batch MAP@10</td><td>0.016 0.018 0.020 0.022 0.024 0.026</td><td>0</td><td>5</td><td colspan="2">10 hours uniform sampling 15 uniform SOGram (α = 0. 01) 20 sqrt sampling sqrt SOGram (α = 0. 01) linear sampling linear SOGram (α = 0. 1)</td><td>Batch MAP@10</td><td>0.06 0.07 0.08 0.09</td><td>0</td><td>5</td><td>10</td><td colspan="2">15 hours uniform sampling 20 25 uniform SOGram (α = 0. 01) 30 sqrt sampling sqrt SOGram (α = 0. 1) linear sampling linear SOGram (α = 0. 001)</td><td>Batch MAP@10</td><td>0.08 0.10 0.12 0.14</td><td>0</td><td>10</td><td>20</td><td>30 uniform sampling 40 hours SOGram uniform (α = 0. 01) sqrt sampling sqrt SOGram (α = 0. 01) linear sampling SOGram linear (α = 0. 001)</td><td>50</td></tr><tr><td colspan="3">language</td><td/><td>uniform sampling</td><td colspan="3">uniform SOGram</td><td/><td/><td colspan="2">sqrt sampling</td><td>sqrt SOGram</td><td/><td/><td/><td>linear sampling</td><td/><td>linear SOGram</td></tr><tr><td colspan="3">simple</td><td/><td>0.0255</td><td colspan="4">0.0266 (+4.2%)</td><td/><td>0.0247</td><td/><td colspan="2">0.0268 (+8.3%)</td><td/><td/><td>0.0272</td><td/><td>0.0300 (+10.7%)</td></tr><tr><td colspan="3">fr</td><td/><td>0.1056</td><td colspan="4">0.1194 (+13.0%)</td><td/><td>0.1047</td><td/><td colspan="2">0.1144 (+9.2%)</td><td/><td/><td>0.1004</td><td/><td>0.1154 (+15.0%)</td></tr><tr><td colspan="3">en</td><td/><td>0.1586</td><td colspan="4">0.1743 (+9.9%)</td><td/><td>0.1543</td><td/><td colspan="3">0.1723 (+11.7%)</td><td/><td>0.1504</td><td/><td>0.1797 (19.5%)</td></tr></table>

Final validation MAP on each dataset, and relative improvement compared to the baselines.

Table content: <table/>

Final validation MAP on Wikipedia en, with batch size |B| = 128.

Table content: <table><tr><td/><td>linear</td><td>linear SOGram</td><td>linear SOGram</td><td>linear SOGram</td></tr><tr><td/><td>sampling</td><td>(α = 0.1)</td><td>(α = 0.01)</td><td>(α = 0.001)</td></tr><tr><td>η = 0.02</td><td>0.0519</td><td colspan="3">0.1331 (+156.4%) 0.1198 (+130.7%) 0.1295 (+149.3%)</td></tr><tr><td>η = 0.01</td><td>0.1371</td><td>0.1329 (-3.0%)</td><td>0.1505 (+9.8%)</td><td>0.1737 (+26.7%)</td></tr><tr><td>η = 0.005</td><td>0.1346</td><td>0.1323 (-1.7%)</td><td>0.1299 (-3.5%)</td><td>0.1569 (+16.6%)</td></tr></table>

Final validation MAP on Wikipedia en, with batch size |B| = 512.

Table content: <table/>

Table 5 gives a basic description of the data size. Note that it is comparable to the simple dataset in the Wikipedia experiments.

Table content: <table><tr><td>Dataset</td><td colspan="3"># users # movies # ratings</td></tr><tr><td>MovieLens</td><td>72K</td><td>10K</td><td>10M</td></tr></table>

Corpus size of the MovieLens dataset.

Table content: <table/>

Peer review (iclr2019):
 This paper proposes a method for estimating non-linear similarities between items using Gramian estimation. This is achieved by having two separate neural networks defined for each item to be compared, which are then combined via a dot product. The proposed innovation in this paper is to use Gramian estimation for the penalty parameter of the optimization which allows for the non-linear case. Two algorithms are proposed which allow for estimation in the stochastic / online setting. Experiments are presented which appear to show good performance on some standard benchmark tasks. Overall, I think this is an interesting set of ideas for an important problem. I have two reservations. First, the organization of the paper needs to be addressed in order to aid user readability. The paper often jumps across sections without giving motivation or connecting language. This will limit the audience of the paper and the work. Second (and more importantly), I found the experiments to be slightly underwhelming. The hyperparameters (batch size, learning rate) and architecture don’t have any rationale attached to them. It is also not entirely clear whether the chosen comparison methods fully constitute the current state of the art. Nonetheless, I think this is an interesting idea and strong work with compelling results. Editorial comments: The organization of this paper leaves something to be desired. The introductions ends very abruptly, and then appears to begin again after the related work section. From what I can tell the first three sections all constitute the introduction and should be merged with appropriate edits to make the narrative clear. “where x and y are nodes in a graph and the similarity is wheter an edge” → typo and sentence ends prematurely.