Year: 2019

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.

heading: 1 INTRODUCTION
text: 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.

heading: 2 PRELIMINARIES

heading: 2.1 NOTATION AND PROBLEM FORMULATION
text: We consider models that consist of two embedding functions u : Rd×X → Rk and v : Rd×Y → Rk, which map a parameter vector1 θ ∈ Rd and feature vectors x, y to embeddings u(θ, x), v(θ, y) ∈ Rk. 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 Rk. 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.3
We denote the training set by T = {(xi, yi, si) ∈ X ×Y ×R}i∈{1,...,n}, where xi, yi are the feature vectors and si is the target similarity for example i. To make notation more compact, we will use ui(θ), vi(θ) as a shorthand for u(θ, xi), v(θ, yi), respectively. As discussed in the introduction, we also assume that we are given a low-similarity prior pij ∈ R for all pairs (i, j) ∈ {1, . . . , n}2. Given a differentiable scalar loss function ` : R× R→ R, the objective function is given by
min θ∈Rd
1
n n∑ i=1 ` (〈ui(θ), vi(θ)〉 , si) + λ n2 n∑ i=1 n∑ j=1 (〈ui(θ), vj(θ)〉 − pij)2, (1)
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 pij 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.
1In 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.
2This also includes cosine similarity models when the embedding functions u, v are normalized. 3We 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).

heading: 2.2 GRAMIAN FORMULATION
text: 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θ ∈ Rn×k the matrix of all left embeddings such that ui(θ) is the i-th row of Uθ, and similarly for Vθ ∈ Rn×k. Then denoting the matrix inner-product by 〈A,B〉 = ∑i,j AijBij , we can rewrite the double sum in (1) as:
g(θ) := 1
n2 n∑ i=1 n∑ j=1 (UθV > θ ) 2 ij = 1 n2 〈 UθV > θ , UθV > θ 〉 . (2)
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
Gu(θ) := 1
n U>θ Uθ =
1
n n∑ i=1 ui(θ)⊗ ui(θ), Gv(θ) := 1 n V >θ Vθ = 1 n n∑ i=1 vi(θ)⊗ vi(θ), (3)
the penalty term becomes g(θ) = 〈Gu(θ), Gv(θ)〉 . (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.

heading: 3 TRAINING EMBEDDING MODELS USING GRAMIAN ESTIMATES
text: 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
fi(θ) := ` (〈ui(θ), vi(θ)〉 , si) (5)
gi(θ) := 1
2 [〈ui(θ), Gv(θ)ui(θ)〉+ 〈vi(θ), Gu(θ)vi(θ)〉], (6)
then (1) is equivalent to minθ∈Rd 1n ∑n i=1[fi(θ) + λgi(θ)]. Intuitively, for each example i, −∇fi(θ) pulls the embeddings ui and vi close to each other (assuming a high similarity si), while −∇gi(θ) creates a repulsive force between ui and all embeddings {vj}j∈{1,...,n}, and between vi and all {uj}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 Ĝu, Ĝv of Gu(θ), Gv(θ). Then the gradient of g(θ) can be approximated by the gradient (w.r.t. θ) of
ĝi(θ, Ĝu, Ĝv) := 〈 ui(θ), Ĝvui(θ) 〉 + 〈 vi(θ), Ĝuvi(θ) 〉 , (7)
as stated in the following proposition.
Proposition 1. If i is drawn uniformly in {1, . . . , n}, and Ĝu, Ĝv are unbiased estimates of Gu(θ), Gv(θ) and independent of i, then ∇θĝi(θ, Ĝu, Ĝv) is an unbiased estimate of∇g(θ).
In a mini-batch setting, one can further average ĝ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 Ĝu, Ĝv , and discuss their tradeoffs. Since each Gramian can be written as a sum of rank-one terms, e.g. Gu(θ) = 1n ∑n i=1 ui(θ)⊗ ui(θ), 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.

heading: 3.1 STOCHASTIC AVERAGE GRAMIAN
text: 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 ûi, v̂i 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 Ŝ(t)u = 1 n ∑n i=1 û (t) i ⊗ û (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
Ĝ(t)u = Ŝ (t) u + β[ui(θ (t))⊗ ui(θ(t))− û(t)i ⊗ û (t) i ], (8)
and similarly for Ĝ(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: β = 1n , 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 Sk+. Proposition 2. Suppose β = 1n in (8). Then for all t, Ĝ (t) u , Ĝ (t) v ∈ Sk+.
Proposition 3. Suppose β = 1 in (8). Then for all t, Ĝ(t)u is an unbiased estimate of Gu(θ(t)).
While taking β = 1 gives an unbiased estimate, note that it does not guarantee that the estimates remain in Sk+. In practice, this can cause numerical issues, but can be avoided by projecting Ĝu, Ĝv on Sk+, using their eigenvalue decompositions. The per-iteration cost of maintaining the Gramian estimates isO(k) to update the caches,O(k2) to update the estimates Ŝu, Ŝv, Ĝu, Ĝv , andO(k3) for projecting on Sk+. Given the small size of the embedding dimension k, O(k3) remains tractable. The memory cost isO(nk), since each embedding needs to be cached (plus a negligibleO(k2) 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 {(xi, yi, si)}i∈{1,...,n}, learning rate η > 0. 2: Initialization phase 3: draw θ randomly 4: ûi ← ui(θ), v̂i ← vi(θ) ∀i ∈ {1, . . . , n} 5: Ŝu ← 1n ∑n i=1 ûi ⊗ ûi, Ŝv ← 1n ∑n i=1 v̂i ⊗ v̂i
6: repeat 7: Update Gramian estimates (i ∼ Uniform(n)) 8: Ĝu ← Ŝu + β[ui(θ)⊗ ui(θ)− ûi ⊗ ûi], Ĝv ← Ŝv + β[vi(θ)⊗ vi(θ)− v̂i ⊗ v̂i] 9: Update model parameters then update caches (i ∼ Uniform(n))

heading: 10: θ ← θ − η∇θ[fi(θ) + λĝi(θ, Ĝu, Ĝv)]

heading: 11: Ŝu ← Ŝu + 1n [ui(θ)⊗ ui(θ)− ûi ⊗ ûi], Ŝv ← Ŝv + 1n [vi(θ)⊗ vi(θ)− v̂i ⊗ v̂i]

heading: 12: ûi ← ui(θ), v̂i ← vi(θ)
text: 13: until stopping criterion

heading: 3.2 STOCHASTIC ONLINE GRAMIAN
text: 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 Ĝu, Ĝv ∈ Sk+, and they seek to minimize the respective losses{
LĜu,Ĝv1 (θ) = 1 n ∑n i=1[fi(θ) + λĝi(θ, Ĝu, Ĝv)] Lθ2(Ĝu, Ĝv) = 1 2‖Ĝu −Gu(θ)‖2F + 12‖Ĝv −Gv(θ)‖2F ,
(9)
where ĝ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. (θ, Ĝu, Ĝv) ∈ Rd × Sk+ × Sk+ is a first-order stationary point for (9) if and only if θ is a first-order stationary point for problem (1) and Ĝu = Gu(θ), Ĝv = Gv(θ).
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∇ĜuL θ 2(Ĝu, Ĝv) = Ĝu−Gu(θ) and can be estimated by Ĝu − ui(θ)⊗ ui(θ), resulting in the update
Ĝ(t)u = (1− α)Ĝ(t−1)u + αui(θ(t))⊗ ui(θ(t)), (10)
and similarly for Ĝ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 Sk+, without the need to project. The per-iteration cost of updating the estimates is O(k2), and the memory cost is O(k2) for storing the Gramians, which are both negligible.
Algorithm 2 SOGram (Stochastic Online Gramian) 1: Input: Training data {(xi, yi, si)}i∈{1,...,n}, learning rates η > 0, α ∈ (0, 1). 2: Initialization phase 3: draw θ randomly, Ĝu, Ĝv ← 0k×k 4: repeat 5: Update Gramian estimates (i ∼ Uniform) 6: Ĝu ← (1− α)Ĝu + αui(θ)⊗ ui(θ), Ĝv ← (1− α)Ĝv + αvi(θ)⊗ vi(θ) 7: Update model parameters (i ∼ Uniform) 8: θ ← θ − η∇θ[fi(θ) + λĝi(θ, Ĝu, Ĝ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, Ĝ(t)u = ∑t τ=1 α(1− α)t−τuiτ (θ(τ))⊗ uiτ (θ(τ)). 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 Ḡ(t)u = ∑t τ=1 α(1−α)t−τGu(θ(τ)). Suppose that there exist σ, δ > 0 such that for all t, Ei ‖ui(θ(t))⊗ ui(θ(t))−Gu(θ(t))‖2F ≤ σ2 and ‖Gu(θ(t+1))−Gu(θ(t))‖F ≤ δ. Then ∀t, E ‖Ĝ(t)u − Ḡ(t)u ‖2F ≤ σ2 α
2− α (11)
‖Ḡ(t)u −G(t)u ‖F ≤ δ(1/α− 1) + (1− α)t‖G(t)u ‖F . (12)
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, Ĝ(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.

heading: 3.3 COMPARISON WITH EXISTING STOCHASTIC METHODS
text: 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
g̃(θ) = 1 |B||B′| ∑ i∈B ∑ j∈B′ µiνj 〈ui(θ), vj(θ)〉2 , (13)
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
g̃(θ) = 〈 1 |B| ∑ i∈B µiui(θ)⊗ ui(θ), 1 |B′| ∑ j∈B′ νjvj(θ)⊗ vj(θ) 〉 . (14)
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 Rk, 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(k2 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.

heading: 4 EXPERIMENTS
text: 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.

heading: 4.1 EXPERIMENTAL SETUP
text: 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 = (xid, xngrams, xcats), where xid is (a one-hot encoding of) the page URL, xngrams is a bag-of-words representation of the set of n-grams of the page’s title, and xcats 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.
4The 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).
Training methods The model is trained using a squared error loss, `(s, s′) = 12 (s− s′)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.

heading: 4.2 QUALITY OF GRAMIAN ESTIMATES
text: 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 Gu(θ(t)), Gv(θ(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 ui(θ(t)), vi(θ(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 ‖Ĝ (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.

heading: 4.3 IMPACT ON TRAINING SPEED AND GENERALIZATION QUALITY
text: 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 summarizes the relative improvement of the final validation MAP.
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.

heading: 5 CONCLUSION
text: 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 Brébisson & Vincent, 2016).

heading: A PROOFS
text: Proof of Proposition 1. Starting from the expression (6) of g(θ) = 〈Gu(θ), Gv(θ)〉, and applying the chain rule, we have
∇g(θ) = ∇〈Gu(θ), Gv(θ)〉 = Ju(θ)[Gv(θ)] + Jv(θ)[Gu(θ)], (15)
where Ju(θ) denotes the Jacobian of Gu(θ), an order-three tensor given by
Ju(θ)l,i,j = ∂Gu(θ)i,j
∂θl , l ∈ {1, . . . , d}, i, j ∈ {1, . . . , n}, and Ju(θ)[Gv(θ)] denotes the vector [ ∑ i,j Ju(θ)l,i,jGv(θ)i,j ]l∈{1,...,d}.
Observing that ĝi(θ, Ĝu, Ĝv) = 〈 Ĝu, ui(θ)⊗ ui(θ) 〉 + 〈 Ĝv, vi(θ)⊗ vi(θ) 〉 , and applying the
chain rule, we have ∇θĝi(θ, Ĝu, Ĝv) = Ju,i(θ)[Ĝv] + Jv,i(θ)[Ĝu], (16) where Ju,i(θ) is the Jacobian of ui(θ)⊗ ui(θ), and
E i∼Uniform
[Ju,i(θ)] = 1
n n∑ i=1 Ju,i(θ) = Ju(θ),
an similarly for Jv,i. We conclude by taking expectations in (16) and using assumption that Ĝu, Ĝv are independent of i.
Proof of Proposition 2. From (8) and the definition of Ŝ(t)u , we have,
Ĝ(t)u = 1
n ∑ j 6=i û (t) j ⊗ û (t) j + 1 n ui(θ (t))⊗ ui(θ(t)),
which is a sum of matrices in the PSD cone Sk+.
Proof of Proposition 3. Denoting by (Ft)t≥0 the filtration generated by the sequence (θ(t))t≥0, and taking conditional expectations in (8), we have
E[Ĝ(t)u |Ft] = Ŝ(t)u + E i∼Uniform [ui(θ (t))⊗ ui(θ(t))− û(t)i ⊗ û (t) i |Ft]
= Ŝ(t)u + 1
n n∑ i=1 [ui(θ (t))⊗ ui(θ(t))− ûi ⊗ ûi]
= 1
n n∑ i=1 ui(θ (t))⊗ ui(θ(t)) = Gu(θ(t)).
Proof of Proposition 4. (θ, Ĝu, Ĝv) ∈ Rd × Sk+ × Sk+ is a first-order stationary point of the game if and only if
∇f(θ) + λ(Ju(θ)[Ĝv] + Jv(θ)[Ĝu]) = 0 (17)〈 Ĝu −Gu(θ), G′ − Ĝu 〉 ≥ 0, ∀G′ ∈ Sk+ (18)〈
Ĝv −Gv(θ), G′ − Ĝv 〉 ≥ 0, ∀G′ ∈ Sk+ (19)
The second and third conditions simply states that ∇ĜuL θ 2(Ĝu, Ĝv) and ∇ĜvL θ 2(Ĝu, Ĝv) define supporting hyperplanes of Sk+ at Ĝu, Ĝv , respectively.
Since Gu(θ) ∈ Sk+, condition (18) is equivalent to Ĝu = Gu(θ) (and similarly, (19) is equivalent to Ĝv = Gv(θ)). Using the expression (15) of ∇g, we get that (17-19) is equivalent to ∇f(θ) + λ∇g(θ) = 0.
Proof of Proposition 5. We start by proving the first bound (11). As stated in Section 3.2, we have, by induction on t, Ĝ(t)u = ∑t τ=1 at−τuiτ (θ
(t))⊗ uiτ (θ(t)), where aτ = α(1− α)τ . And by definition of Ḡ(t), we have Ḡ(t)u = ∑t τ=1 at−τGu(θ (τ)). Thus,
Ĝ(t)u − Ḡ(t)u = t∑
τ=1
at−τ∆ (τ) u
where ∆(τ)u = uiτ (θ (τ))⊗ uiτ (θ(τ))−Gu(θ(τ)) are zero-mean random variables. Thus, taking the second moment, and using the first assumption (which simply states that the variance of ∆(τ)u is bounded by σ2), we have
E ‖Ĝ(t)u − Ḡ(t)u ‖2F = E ∥∥∥∥∥ t∑ τ=1 at−τ∆ (τ) u ∥∥∥∥∥ 2
F
= t∑ τ=1 a2t−τ E ‖∆(τ)u ‖2F
≤ σ2α2 t−1∑ τ=0 (1− α)2τ = σ2α2 1− (1− α) 2t 1− (1− α)2 ≤ σ2 α 2− α,
which proves the first inequality (11).
To prove the second inequality, we start from the definition of Ḡ(t)u :
‖Ḡ(t)u −G(t)u ‖F = ‖ t∑
τ=1
at−τ (G (τ) u −G(t)u )− (1− α)tG(t)u ‖F
≤ t∑
τ=1
at−τ‖G(τ)u −G(t)u ‖F + (1− α)t‖G(t)u ‖F , (20)
where the first equality uses that fact that ∑t τ=1 at−τ = 1 − (1 − α)t. Focusing on the first term, and bounding ‖G(τ)u −G(t)u ‖F ≤ (t− τ)δ by the triangle inequality, we get t∑
τ=1
at−τ‖G(τ)u −G(t)u )‖F ≤ δ t∑
τ=1 at−τ (t− τ) = δα t−1∑ τ=0 τ(1− α)τ
= δα(1− α) d dα [ − t−1∑ τ=0 (1− α)τ ]
= δα(1− α) d dα
[ −1− (1− α) t
α ] ≤ δα(1− α) 1
α2 . (21)
Combining (20) and (21), we get the desired inequality (12).

heading: B GENERALIZATION TO LOW-RANK PRIORS
text: So far, we have assumed a uniform zero prior to simplify the notation. In this section, we relax this assumption. Suppose that the prior is given by a low-rank matrix P = QR>, where Q,R ∈ Rn×kP . In other words, the prior for a given pair (i, j) is given by the dot product of two vectors pij = 〈qi, rj〉. In practice, such a low-rank prior can be obtained, for example, by first training a simple low-rank matrix approximation of the target similarity matrix.
Given this low-rank prior, the penalty term (2) becomes
gP (θ) = 1
n2 n∑ i=1 n∑ j=1 [UθV > θ −QR>]2ij
= 1 n2 〈 UθV > θ −QR>, UθV >θ −QR> 〉 = 1
n2 [〈 U>θ Uθ, V > θ Vθ 〉 − 2 〈 U>θ Q,V > θ R 〉 + c ] ,
where c = 〈 Q>Q,R>R 〉 is a constant that does not depend on θ. Here, we used a superscript P in
gP to disambiguate the zero-prior case.
Now, if we define weighted embedding matrices{ Hu(θ) := 1 nUθQ = 1 n ∑n i=1 ui(θ)⊗ qi
Hv(θ) := 1 nVθR = 1 n ∑n i=1 vi(θ)⊗ ri,
the penalty term becomes
gP (θ) = 〈Gu(θ), Gv(θ)〉 − 2 〈Hu(θ), Hv(θ)〉+ c.
Finally, if we maintain estimates Ĥu, Ĥv of Hu(θ), Hv(θ), respectively (using the methods proposed in Section 3), we can approximate∇gP (θ) by the gradient of
ĝPi (θ, Ĝu, Ĝv, Ĥu, Ĥv) :=〈 ui(θ), Ĝvui(θ) 〉 + 〈 vi(θ), Ĝuvi(θ) 〉 − 2 〈 ui(θ), Ĥvqi 〉 − 2 〈 vi(θ), Ĥuri 〉 . (22)
Proposition 1 and Algorithms 1 and 2 can be generalized to the low-rank prior case by adding updates for Ĥu, Ĥv , and by using expression (22) of ĝPi when computing the gradient estimate.
Proposition 6. If i is drawn uniformly in {1, . . . , n}, and Ĝu, Ĝv , Ĥu, Ĥv are unbiased estimates of Gu(θ), Gv(θ), Hu(θ), Hv(θ), respectively, then ∇θĝPi (θ, Ĝu, Ĝv, Ĥu, Ĥv) is an unbiased estimate of∇gP (θ).
Proof. Similar to the proof of Proposition 1.
The generalized versions of SAGram and SOGram are stated below, where the differences compared to the zero prior case are highlighted. Note that, unlike the Gramian matrices, the weighted embedding matrices Hu, Hv are not symmetric, thus we do not project their estimates.
Algorithm 3 SAGram (Stochastic Average Gramian) with low-rank prior 1: Input: Training data {(xi, yi, si)}i∈{1,...,n}, low-rank priors {qi, ri}i∈{1,...,n} 2: Initialization phase 3: draw θ randomly 4: ûi ← ui(θ), v̂i ← vi(θ) ∀i ∈ {1, . . . , n} 5: Ŝu ← 1n ∑n i=1 ûi ⊗ ûi, Ŝv ← 1n ∑n i=1 v̂i ⊗ v̂i
6: T̂u ← 1n ∑n i=1 ûi ⊗ qi, T̂v ← 1n ∑n
i=1 v̂i ⊗ ri 7: repeat 8: Update Gramian estimates (i ∼ Uniform(n)) 9: Ĝu ← Ŝu + β[ui(θ)⊗ ui(θ)− ûi ⊗ ûi], Ĝv ← Ŝv + β[vi(θ)⊗ vi(θ)− v̂i ⊗ v̂i] 10: Update weighted embedding estimates 11: Ĥu ← T̂u + λ[(ui(θ)− ûi)⊗ qi] 12: Ĥv ← T̂v + λ[(vi(θ)− v̂i)⊗ ri] 13: Update model parameters then update caches (i ∼ Uniform(n))

heading: 14: θ ← θ − η∇θ[fi(θ) + λp̃ĝP (θ, Ĝu, Ĝv, Ĥu, Ĥv)]

heading: 15: Ŝu ← Ŝu + 1n [ui(θ)⊗ ui(θ)− ûi ⊗ ûi], Ŝv ← Ŝv + 1n [vi(θ)⊗ vi(θ)− v̂i ⊗ v̂i]
text: 16: T̂u ← T̂u + 1n [(ui(θ)− ûi)⊗ qi], T̂v ← T̂v + 1n [(vi(θ)− v̂i)⊗ ri] 17: ûi ← ui(θ), v̂i ← vi(θ) 18: until stopping criterion
Algorithm 4 SOGram (Stochastic Online Gramian) with low-rank prior 1: Input: Training data {(xi, yi, si)}i∈{1,...,n}, low-rank priors {qi, ri}i∈{1,...,n} 2: Initialization phase 3: draw θ randomly 4: Ĝu, Ĝv ← 0k×k 5: repeat 6: Update Gramian estimates (i ∼ Uniform(n)) 7: Ĝu ← (1− α)Ĝu + αui(θ)⊗ ui(θ), Ĝv ← (1− α)Ĝv + αvi(θ)⊗ vi(θ) 8: Update weighted embedding estimates 9: Ĥu ← (1− α)Ĥu + αui(θ)⊗ qi, Ĥv ← (1− α)Ĥv + αvi(θ)⊗ ri 10: Update model parameters (i ∼ Uniform(n))

heading: 11: θ ← θ − η∇θ[fi(θ) + λp̃ĝP (θ, Ĝu, Ĝv, Ĥu, Ĥv)]
text: 12: until stopping criterion

heading: C NON-UNIFORM WEIGHTS
text: 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
gW (θ) = 1
n2 n∑ i=1 n∑ j=1 aibj 〈ui(θ), vj(θ)〉2 ,
where ai, bj are positive left and right weights, respectively. Here we used a superscript W in gW to disambiguate the uniform-weight case. Then using a similar transformation to Section 2.2, we can rewrite gW as follows:
gW (θ) =
〈 1
n n∑ i=1 aiui ⊗ ui, 1 n n∑ j=1
bjvj ⊗ vj 〉 ,
i.e. gW 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 ui ⊗ ui by the appropriate embedding weight ai (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 n2 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
g(θ) = 〈Gu(θ), Gv(θ)〉 = 〈 1
n n∑ i=1
ui(θ)⊗ ui(θ), Gv(θ) 〉 = 1
n n∑ i=1 〈ui(θ), Gv(θ)ui(θ)〉 , (23)
which is a quadratic form in the left embeddings ui (and similarly for vj , by symmetry). In particular, the partial derivative of the Gramian term with respect to an embedding ui is
∂g(θ)
∂ui =
2 n Gv(θ)ui(θ) = 2 n  1 n n∑ j=1 vj(θ)⊗ vj(θ) ui(θ). Each term (vj ⊗ vj)ui = vj 〈vj , ui〉 is simply the projection of ui on vj (scaled by ‖vj‖2). Thus the gradient of g(θ) with respect to ui is an average of scaled projections of ui on each of the right embeddings vj , and moving in the direction of the negative gradient simply moves ui 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 〈 ui(θ (t)), vj(θ (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 (λ = 103), 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(gu(θ), 1/k)+max(gv(θ), 1/k), where k is the embedding dimension, gu(θ) = 1n2 ∑n i=1 ∑n j=1 〈ui, uj〉
2, and similarly for gv . They optimize this term using candidate sampling. We can also apply the same Gramian transformation as in Section 2.2 to write gu(θ) = 〈Gu(θ), Gu(θ)〉, and gv(θ) = 〈Gv(θ), Gv(θ)〉, 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.

heading: E FURTHER EXPERIMENTS ON WIKIPEDIA
text: 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(θ)
∂ui =
2 n GV (θ)ui(θ)
∂ĝ(θ)
∂ui =
2 n ĜV ui(θ),
the estimation error of the gradient with respect to the left embeddings u is
‖∇ug −∇uĝ‖22 = 4
n2 n∑ i=1 ‖(Gv − Ĝv)ui‖2
= 4
n2 n∑ i=1 〈 (Gv − Ĝv)ui, (Gv − Ĝv)ui 〉 = 4
n2 n∑ i=1 〈 Gv − Ĝv, (Gv − ĜV )ui ⊗ ui 〉 = 4
n
〈 Gv − Ĝv, (Gv − Ĝv)Gu 〉 .
This last expression can be interpreted as a Frobenius norm of the right Gramian estimation error Gv − Ĝv , weighted by the left Gramian Gu, thus the gradient error is closely related to the Gramian error. Figure 6 shows the gradient estimation quality on Wikipedia simple, measured by the normalized squared norm ‖∇uĝ−∇ug‖2‖∇ug‖2 . The results are similar to the Gramian estimation errors reported in Figure 2.
E.2 EFFECT OF THE SAMPLING DISTRIBUTION
Comparing the different baselines methods on simple (Figure 3), we observe that uniform sampling performs better than sqrt, despite having a worse Gramian estimate according to Figure 2a. One possible explanation is that the sampling distribution affects both the quality of the Gramian estimates, and the frequency at which the item embeddings are updated, which in turn affects the MAP. In particular, tail items are updated more frequently under uniform than other distributions, and this may have a positive impact on the MAP.
E.3 EFFECT OF GRAMIAN LEARNING RATE α AND BIAS-VARIANCE TRADEOFF
In addition to the experiments of Section 4, we also evaluated the effect of the Gramian learning rate α on the quality of the Gramian esimates and generalization performance on Wikipedia en. Figure 7 shows the validation MAP of the SOGram method for different values of α (together with the basline for reference). This reflects the bias-variance tradeoff dicussed in Proposition 5: with a lower α, progress is initially slower (due to the bias introduced in the Gramian estimates), but the final performance is better. Given a limited training time budget, this suggests that a higher α can be preferable.
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 Ĝu and (the large sample approximation of) the true Gramian Gu. 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.

heading: F EXPERIMENT ON MOVIELENS DATA
text: 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. Table 5 gives a basic description of the data size. Note that it is comparable to the simple dataset in the Wikipedia experiments.
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′) = 12 (s − s′)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
estimation of g(θ) induces better regularization. The impact on training speed is also remarkable in this case, SOGram with α = 0.1 achieves a better validation performance in under 1 hour of training than the sampling baseline in 6 hours.

Review summaries:
This paper proposes an efficient algorithm to learn  neural embedding models with a dot-product structure over very large corpora. The main method is to reformulate the objective function in terms of generalized Gramiam matrices, and maintain estimates of those matrices in the training process. The algorithm uses less time and achieves significantly better quality than sampling based methods. 

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.