Year: 2019

Abstract:
The leading approaches to tensor completion and link prediction are based on the canonical polyadic (CP) decomposition of tensors. While these approaches were originally motivated by low rank approximations, the best performances are usually obtained for ranks as high as permitted by computation constraints. For large scale factorization problems where the factor dimensions have to be kept small, the performances of these approaches tend to drop drastically. The other main tensor factorization model, Tucker decomposition, is more flexible than CP for fixed factor dimensions, so we expect Tucker-based approaches to yield better performance under strong constraints on the number of parameters. However, as we show in this paper through experiments on standard benchmarks of link prediction in knowledge bases, ComplEx, (Trouillon et al., 2016), a variant of CP, achieves similar performances to recent approaches based on Tucker decomposition on all operating points in terms of number of parameters. In a control experiment, we show that one problem in the practical application of Tucker decomposition to large-scale tensor completion comes from the adaptive optimization algorithms based on diagonal rescaling, such as Adagrad. We present a new algorithm for a constrained version of Tucker which implicitly applies Adagrad to a CP-based model with an additional projection of the embeddings onto a fixed lower dimensional subspace. The resulting Tucker-style extension of ComplEx obtains similar best performances as ComplEx, with substantial gains on some datasets under constraints on the number of parameters.

heading: 1 INTRODUCTION
text: The problems of representation learning and link prediction in multi-relational data can be formulated as a binary tensor completion problem, where the tensor is obtained by stacking the adjacency matrices of every relations between entities. This tensor can then be intrepreted as a "knowledge base", and contains triples (subject, predicate, object) representing facts about the world. Link prediction in knowledge bases aims at automatically discovering missing facts (Bordes et al., 2011; Nickel et al., 2011; Bordes et al., 2013; Nickel et al., 2016a; Nguyen, 2017).
State of the art methods use the canonical polyadic (CP) decomposition of tensors (Hitchcock, 1927) or variants of it (Trouillon et al., 2016; Kazemi & Poole, 2018; Lacroix et al., 2018). While initially motivated by low-rank assumptions on the underlying ground-truth tensor, the best performances are obtained by setting the rank as high as permitted by computational constraints, using tensor norms for regularization (Lacroix et al., 2018). However, for large scale data where computational or memory constraints require ranks to be low (Lerer et al., 2019), performances drop drastically.
Tucker decomposition is another multilinear model which allows richer interactions between entities and predicate vectors. A special case of Tucker decomposition is RESCAL (Nickel et al., 2011), in which the relations are represented by matrices and entities factors are shared for subjects and objects. However, an evaluation of this model in Nickel et al. (2016b) shows that RESCAL lags behind other methods on several benchmarks of interest. Recent work have obtained more competitive results with similar models (Balažević et al., 2019b; Wang et al., 2019), using different regularizers or deep learning heuristics such as dropout and label smoothing. Despite these recent efforts, learning Tucker decompositions remains mostly unresolved. Wang et al. (2019) does not achieve state of the art results on standard benchmarks, and we show (see Figure 3) that the performances reported by
Balažević et al. (2019b) are actually matched by ComplEx (Trouillon et al., 2016; Lacroix et al., 2018) optimized with Adam, which has less hyperparameters.
In this work, we overcome some of the difficulties associated with learning a Tucker model for knowledge base completion. Balažević et al. (2019b) use deep-learning mechanisms such as batch normalization (Ioffe & Szegedy, 2015), dropout (Srivastava et al., 2014) or learning-rate annealing to address both regularization and optimization issues. Our approach is different: We factorize the core tensor of the Tucker decomposition with CP to obtain a formulation which is closer to CP and better understand what difficulties appear. This yields a simple approach, which has a single regularization hyperparameter to tune for a fixed model specification.
The main novelty of our approach is a more careful application of adaptive gradient techniques. State-of-the-art methods for tensor completion use optimization algorithms with adaptive diagonal rescaling such as Adagrad (Duchi et al., 2011) or Adam (Kingma & Ba, 2014). Through control experiments in which our model is equivalent to CP up to a fixed rotation of the embeddings, we show that one of the difficulties in training Tucker-style decompositions can be attributed to the lack of invariance to rotation of the diagonal rescaling. Focusing on Adagrad, we propose a different update rule that is equivalent to implicitely applying Adagrad to a CP model with a projection of the embedding to a lower dimensional subspace.
Combining the Tucker formulation and the implicit Adagrad update, we obtain performances that match state-of-the-art methods on the standard benchmarks and achieve significanly better results for small embedding sizes on several datasets. Compared to the best current algorithm for Tucker decomposition of Balažević et al. (2019b), our approach has less hyperparameters, and we effectively report better performances than the implementation of ComplEx of Lacroix et al. (2018) in the regime of small embedding dimension.
We discuss the related work in the next section. In Section 3, we present a variant of the Tucker decomposition which allows to interpolate between Tucker and CP. The extreme case of this variant, which is equivalent to CP up to a fixed rotation of the embedding, serves as control model to highlight the deficiency of the diagonal rescaling of Adagrad for Tucker-style decompositions in experiments reported in Section 4 . We present the modified version of Adagrad in Section 5 and present experimental results on standard benchmarks of knowledge base completion in Section 7.

heading: 2 LINK PREDICTION IN KNOWLEDGE BASES
text: Notation Tensors and matrices are denoted by uppercase letters. For a matrix U , ui is the vector corresponding to the i-th row of U . The tensor product is written ⊗ and the Hadamard product (i.e., elementwise product) is written ⊙.

heading: 2.1 LEARNING SETUP
text: A knowledge base consists of a set S of triples (subject, predicate, object) that represent (true) known facts. The goal of link prediction is to recover facts that are true but not in the database. The data is represented as a tensor X̃ ∈ {0, 1}N×L×N for N the number of entities and L the number of predicates. Given a training set of triples, the goal is to provide a ranking of entities for queries of the type (subject, predicate, ?) and (?, predicate, object). Following Lacroix et al. (2018), we use the cross-entropy as a surrogate of the ranking loss. As proposed by Lacroix et al. (2018) and Kazemi & Poole (2018), we include reciprocal predicates: for each predicate P in the original dataset, and given an item o, each query of the form (?, P , o) is reformulated as a query (o, P−1, ?), where o is now the subject of P−1. This doubles the effective number of predicates but reduces the problem to queries of the type (subject, predicate, ?) only.
For a given triple (i, j, k) ∈ S, the training loss function for a tensor X is then
ℓi,j,k(X) = −Xi,j,k + log ( ∑
k′ ̸=k
exp(Xi,j,k′) ) . (1)
For a tensor decomposition model X(θ) parameterized by θ, the parameters θ̂ are found by minimizing the regularized empirical risk with regularizer Λ:
θ̂ = argmin θ L(θ) = argmin θ
1 |S| ∑
(i,j,k)∈S
ℓi,j,k(X(θ)) + νΛ(θ). (2)
This work studies specific models for X(θ), inspired by CP and Tucker decomposition. We discuss the related work on tensor decompositions and link prediction in knowledge bases below.

heading: 2.2 RELATED WORK

heading: 2.2.1 CANONICAL DECOMPOSITION AND ITS VARIANTS
text: The canonical polyadic (CP) decomposition of a tensor X is defined entrywise by
∀i, j, k, Xi,j,k = ⟨ui, vj , wk⟩ := d∑
r=1
uirvjrwkr.
The smallest value of d for which this decomposition exists is the rank of X . Each element Xi,j,k is thus represented as a multi-linear product of the 3 embeddings in Rd associated respectively to the ith subject, the jth predictate and the kth object.
CP currently achieves near state-of-the-art performances on standard benchmarks of knowledge base completion (Kazemi & Poole, 2018; Lacroix et al., 2018). Nonetheless, the best reported results are with the ComplEx model (Trouillon et al., 2016), which learns complex-valued embeddings and sets the embeddings of the objects to be the complex conjugate of the embeddings of subjects, i.e., wk = ūk. Prior to ComplEx, Dismult was proposed (Yang et al., 2014) as a variant of CP with wk = uk. While this model obtained good performances (Kadlec et al., 2017), it can only model symmetric relations and does not perform as well as ComplEx. CP-based models are optimized with vanilla Adam or Adagrad and a single regularization parameter (Trouillon et al., 2016; Kadlec et al., 2017; Lacroix et al., 2018) and do not require additional heuristics for training.

heading: 2.2.2 TUCKER DECOMPOSITION AND ITS VARIANTS
text: Given a tensor X of size N×L×N , the Tucker decomposition of X is defined entrywise by
∀i, j, k, Xi,j,k = ⟨ui ⊗ vj ⊗ wk, C⟩ := d1∑
r1=1 d2∑ r2=1 d3∑ r3=1 Cr1,r2,r3uir1vjr2wkr3 .
The triple (d1, d2, d3) are the rank parameters of the decomposition. We also use a multilinear product notation X = [[C;U, V,W ]], where U, V,W are the matrices whose rows are respectively uj , vk, wl and C the three dimensional d1 × d2 × d3 core tensor. Note that the CP decomposition is a Tucker decomposition in which d1 = d2 = d3 = d and C is the identity, which we write [[U, V,W ]]. With a non-trivial core tensor, Tucker decomposition is thus more flexible than CP for fixed embedding size. In knowledge base applications, we typically have d ≤ L ≪ N , so the vast majority of the model parameters are in the embedding matrices of the entities U and W . When constraints on the number of model parameters arise (e.g., memory constraints), Tucker models appear as natural candidates to increase the expressivity of the decomposition compared to CP with limited impact on the total number of parameters.
While many variants of the Tucker decomposition have been proposed in the literature on tensor factorization (see e.g., Kolda & Bader, 2009), the first approach based on Tucker for link prediction in knowledge bases is RESCAL (Nickel et al., 2011). RESCAL uses a special form of Tucker decomposition in which the object and subject embeddings are shared, i.e., U = W , and it does not compress the relation matrices. In the multilinear product notation above, a RESCAL model is thus written as X = [[C;U, I, U ]]. Despite some success on a few smaller datasets, RESCAL performances drop on larger datasets (Nickel et al., 2016b). This decrease in performances has been attributed either to improper regularization (Nickel et al., 2011) or optimization issues (Xue et al., 2018). Balažević et al. (2019b) revisits Tucker decomposition in the context of large-scale knowledge bases and resolves some of the optimization and regularization issues using learning rate
annealing, batch-normalization and dropout. It comes at the price of more hyperparameters to tune for each dataset (label smoothing, three different dropouts and a learning rate decay), and as we discuss in our experiments, the results they report are not better than ComplEx for the same number of parameters.
Two methods were previously proposed to interpolate between the expressivity of RESCAL and CP. Xue et al. (2018) expands the HolE model (Nickel et al., 2016b) (and thus the ComplEx model (Hayashi & Shimbo, 2017)) based on cross-correlation of embeddings to close the gap in expressivity with the Tucker decomposition for a fixed embedding size. Jenatton et al. (2012) express the relation matrices in RESCAL as low-rank combination of a family of matrices. We describe the link between these approaches and ours in Appendix 9.4. None of these approach however studied the effect of their formulation on optimization, and reported results inferior to ours.

heading: 2.2.3 OTHER APPROACHES
text: (Graph) neural networks for link prediction Several methods have introduced models that go beyond the form of Tucker and canonical decompositions. ConvE (Dettmers et al., 2018) uses a convolution on a 2D tiling of the subject and relation embeddings as input to a 2-layer neural net that produces a new embedding for the pair, then compares to the object embedding. Graph neural networks (Scarselli et al., 2009; Niepert et al., 2016; Li et al., 2016; Bruna et al., 2014) have recently gained popularity and have been applied to link prediction in knowledge bases by Schlichtkrull et al. (2018). This model uses a graph convolutional architecture to generate a variant of CP.
Poincaré embeddings Poincaré embeddings have been proposed as an alternative to usual tensor decomposition approaches to learn smaller embeddings when the relations are hierarchical (Nickel & Kiela, 2017). The method has recently been extended to link prediction in relational data with very good performance trade-offs for small dimensional embeddings on the benchmark using WordNet (Balažević et al., 2019a), which contains relationships such as hypernyms and hyponyms which are purely hierarchical. However, such good results do not extend to other benchmarks.

heading: 3 INTERPOLATING BETWEEN CP AND TUCKER
text: In order to better understand the underlying difficulties in learning (variants of) Tucker decompositions compared to CP, our analysis starts from a Tucker model in which the core tensor is itself decomposed with CP. Given a N × L × N tensor, a fixed d and assuming a (d, d, d) Tucker decomposition to simplify notation, a Tucker model where the core tensor is itself decomposed with a rank-D CP can be written as (details are given in Appendix 9.3):
Xijk = ⟨ui ⊗ vj ⊗ wk, C⟩ = ⟨P1ui, P2vj , P3wk⟩ or equivalently X = [[UP⊤1 , V P⊤2 ,WP⊤3 ]],
where P1, P2, P3 are all D × d matrices. Since most knowledge bases have much fewer predicates than entities (L≪ N ), the dimension of the predictate factors has little impact on the overall number of model parameters. So in the remainder of the paper, we always consider P2 = I . Learning matrices U, V,W,P1, P3 of this decomposition simultaneously leads to the following model, which we call CP-Tucker (CPT):
(CPT)Xijk = ⟨P1ui, vj , P3wk⟩, ui, wk ∈ Rd, vj ∈ RD, Pi ∈ RD×d.
The CPT model is similar to a CP model except that the embedding matrices U and W have an additional low-rank constraint (d instead of D). We say that the model interpolates between CP and Tucker because for D = d it is equivalent to CP (as long as P1 and P3 are full rank), whereas for D = d2 we recover a full Tucker model because the matrices P1 and P3 can be chosen such that ⟨P1ui, vj , P3wk⟩ = uiMat(vj)vTj , where Mat is the operator that maps a d2 vector to a d×d matrix (see Appendix 9.5).
CPT is similar to CANDELINC (Carroll et al., 1980), except that in CANDELINC the factors U , V and W are fixed and used to compress the data in order to efficiently learn the Pi. Closer to CPT, Bro & Andersson (1998) first learn a Tucker3 decomposition of X before applying CANDELINC using the learned factors. These methods are only applicable to least-square estimation, and for tensors of smaller scale than knowledge bases.
Fixed projection matrices: The Projected Canonical Polyadic (PCP) Decomposition In order to clarify the difficulty that arise when learning a CPT model compared to a CP model, we study a simpler model in which the matrices P1 and P3 are not learned but rather fixed during training and taken as random matrices with orthonormal columns. We call the resulting model the Projected Canonical Polyadic (PCP) decomposition, since P1, P3 project the embeddings of dimension d into a higher dimension D:
(PCP)Xijk = ⟨P1ui, vj , P3wk⟩, ui, wk ∈ Rd, vj ∈ RD, fixed P1, P3 ∈ RD×d
When D = d the matrices Pi are then fixed unitary transformations. The PCP (or CPT) model in this case D = d is then equivalent to a CP model, up to a fixed invertible transformation of the embeddings. The capacity of the model grows beyond that of CP as D increases up to d2.

heading: 4 MOTIVATION: OPTIMIZATION ISSUES WITH CPT AND PCP
text: As discussed in the related works, previous results suggest that Tucker models are more difficult to train than CP models. The goal of this section is to isolate an issue faced with CPT/PCP models when trained with vanilla adaptive gradient methods such as Adagrad or Adam.

heading: 4.1 CONTROL EXPERIMENT: UNITARY P1 AND P3 IN PCP
text: When D = d in PCP, the model becomes equivalent to CP. Indeed, the matrices P1 and P3 are unitary (P1P⊤1 = P3P ⊤ 3 = I) and so [[(UP1)P ⊤ 1 , V, (WP3)P ⊤ 3 ]] = [[U, V,W ]]. There is no practical interest in considering this degenerate case of PCP, we only use it in the following toy experiment to exhibit one of the difficulties encountered when training PCP.
We perform a simple control experiment in which we take one of the standard benchmarks of link prediction in knowledge bases, called FB15K-237, and train a CP model for different values of the rank D and a PCP model with D = d with vanilla Adagrad. The full experimental protocol, including hyperparameter tuning, is similar to our main experiments and is described in Section 7.2. Figure 1a plots the performances in terms of the standard metric mean reciprocal rank (higher is better) as a function of D of CP (blue curve) and PCP (red curve, called PCP (Adagrad)).
We observe that CP obtains significantly better performances than CPT for larger embedding dimension D. Since in this toy experiment CP and PCP can represent exactly the same tensors and have equivalent regularizations, the only difference between the algorithms that can explain the difference in performances is in how the optimization is carried out, namely the diagonal rescaling performed by Adagrad: Adagrad adapts the learning rate on a per-parameter basis, depending on previous and
current gradients, and is therefore not invariant by the addition of the matrices P1 and P2 even if these are unitary (we provide the formal justification in the next section). This is shown experimentally in Figure 1b where we plot the average Adagrad coefficients for each embedding dimensions (i.e., adagrad coefficients of subject/object embedding matrices averaged by column). The addition of the random P1 and P2 flattens the Adagrad weights, which in turn removes all the benefit of the adaptive rescaling of the algorithm.
For reference, we also tried to directly learn all parameters including P1 and P3 (i.e., learn a CPT model) with vanilla Adagrad. The performances obtained are also lower than those of CP, as shown in Figure 1a (orange curve).
5 A ROTATION INVARIANT ADAGRAD: ADAimp
In this section, we study the optimization problem in more details, and more precisely the effect of the diagonal rescaling performed by Adagrad. As a remainder, given a sequence of stochastic gradients g(t) of L and denoting G(t) = ϵI + ∑t τ=1 g (τ)g(τ) ⊤ , the (practical) AdaGrad update is:
θ(t+1)p = θ (t) p − ηg(t)p /
√ G (t) pp or equivalently θ(t+1) = θ(t) − η Diag(G(t))−1/2
where Diag(G) is the diagonal matrix obtained by extracting the diagonal elements of G.

heading: 5.1 TWO EQUIVALENT PARAMETRIZATIONS OF PCP
text: The following decomposition is equivalent to PCP, but its embeddings are expressed in RD:
(PCPFULL)Xi,j,k = ⟨P1P⊤1 ui, vj , P3P⊤3 wk⟩, with ui, vj , wk ∈ RD, fixed P1, P2 ∈ RD×d.
Note that with ui = P ⊤ 1 ui and wk = P ⊤ 3 wk, we go from PCPfull to PCP. The practical differences between PCP and PCPfull are that PCPfull learn embeddings in the high-dimensional space, maintaining the low-rank structure of the overall entity embeddings through the orthogonal projections P1P T 1 and P3P T 3 . The practical interest of PCPfull is not in terms of modeling but rather from the optimization perspective with Adagrad because it has a structure that is closer to that of CP.
Indeed, for d = D, P1 and P3 disappear in PCPfull, so that optimizing a PCPfull model with Adagrad is equivalent to optimizing a CP model with Adagrad. This property suggests an alternative to the vanilla PCP + Adagrad algorithm, which we call implicit Adagrad:
Implicit Adagrad: Adaimp The approach we propose is to effectively optimize PCPfull with Adagrad. However, when d < D, which is the interesting case for PCP, we notice that we do not need to maintain embeddings in RD. Our approach, called Adaimp, computes the gradients and Adagrad coefficients with respect to ui, wk ∈ RD, but the full dimension factor matrices U and W are never explicitly stored in memory. Rather, we store ui = P ⊤ 1 ui and wk = P ⊤ 3 wk ∈ Rd, which is all that is required for any model computation in PCPfull since P1 and P3 are fixed. Overall, the effective model parameters are exactly the same as in PCP, and we call this approach PCP +Adaimp.
An Adaimp update is described in Algorithm 4. While PCP +Adagrad and PCP +Adaimp work with the same number of model parameters, the fundamental difference is the computation of Adagrad coefficients. Since Adaimp effectively applies Adagrad to PCPfull, we need to maintain the Adagrad coefficients in RD even when d < D: the overall update is first computed in RD and projected back to Rd after the application of Adagrad rescaling. In constrat, in vanilla PCP +Adagrad, the gradient is projected to Rd before Adagrad rescaling.
5.2 IMPLICIT OPTIMIZATION OF PCPFULL : ADAimp
In this section, we discuss more formally the Adaimp updates, and how they compare to PCP +Adagrad. In the following, Ũ , W̃ are in RN×d, whereas U, V and W are in RN×D. Using d ≤ D, P1, P3 ∈ RD×d, and we use the notation Π1 = P1P⊤1 and Π3 = P3P⊤3 .
The empirical risk L can be expressed as a function of three matrices M (1), M (2) and M (3) corresponding to factors of a CP decomposition. We focus on the following problems :
(PCP) argmin Ũ,V,W̃ L(ŨP⊤1 , V, W̃P⊤3 ) (PCPFULL) argmin U,V,W L(UΠ⊤1 , V,WΠ⊤3 )
We focus on a step at time (t) on vectors ũi and ui. We assume that at this time t, the tensor iterates are the same, that is Ũ = UP⊤1 (resp. W̃ = WP ⊤ 1 ) so that [[ŨP ⊤ 1 , V, W̃P ⊤ 3 ]] = [[UΠ ⊤ 1 , V,WΠ ⊤ 3 ]]. In this case, the gradient ∇ M (1) i L is the same in both optimization problems, we denote it by g(t)i .
Let G(t)i = ϵId + ∑t τ=1 g (τ) i g (τ) i ⊤ . The updates for (PCP) are:
ũt+1i = ũ t i − ηDiag(P⊤1 G (t) i P1) −1/2P⊤1 g (t) i . (3)
Note that due to the presence of P1 inside the Diag operator, the update (3) is not rotation invariant. Moreover, for random P1, the matrix P⊤1 G (t) i P1 will be far from diagonal with high probability, making its diagonal meaningless. This is visualized in Figure 1b.
Similar updates can be derived for (PCPfull):
ut+1i = u t i − ηDiag(Π⊤1 G (t) i Π1) −1/2Π⊤1 g (t) i . (4)
As a sanity check, clearly, for d = D and Π1 = I , the update (4) is equivalent to the Adagrad update for the CP model. In the general case d ≤ D, in order to avoid storing U ∈ RN×D, we apply these updates implicitly with Adaimp, by storing u(t)i = P ⊤ 1 u (t) i in Rd. Let us compare the updates :
(ADAGRAD) ũt+1i = ũ t i − η × Adagrad in Rd︷ ︸︸ ︷ Diag(P⊤1 G (t) i P1) −1/2 × projection to Rd︷ ︸︸ ︷ P⊤1 g (t) i
(ADAimp) ut+1i = u t i − η × P⊤1︸︷︷︸
projection to Rd
× Diag(Π⊤1 G (t) i Π1) −1/2︸ ︷︷ ︸ Adagrad in RD
× Π⊤1 g (t) i︸ ︷︷ ︸
projection to Im(Π)∈RD
Going back to our control experiment, we note on Figure 1a that PCP +Adaimp matches the performances of CP+Adagrad for all D, indicating that we fixed this optimization issue.
5.3 ALTERNATIVES TO ADAimp
Another solution would be to use Adagrad projected on the column space of Π, but we show in Appendix 9.1 that even with the diagonal approximation, this is impractical. Note that the version of Adagrad which uses the full matrix G(t) is rotation invariant (see Appendix 9.2 for details), but it cannot be used at the scale of our problems.
It could be argued that the strength of the AdaGrad algorithm in our context mostly comes from its adaptation to the different frequencies of updates of each embedding. In fact, this is one of the examples chosen in Duchi et al. (2011) to display the advantages of AdaGrad compared to stochastic gradient descent. A version of AdaGrad that would only keep one coefficient per embedding (we call this version Adarow) would be invariant to unitary transforms by design and would adapt to the various update frequencies of different embeddings. In fact, this version of AdaGrad is used to save memory in Lerer et al. (2019). We test this algorithm in Appendix 9.8. The difference in performances shows that this adaptation is not sufficient to recover the performances of the finer diagonal AdaGrad in our setting. The claim that the approximations made in Adaimp are indeed better is further backed by the experiments in the next section .

heading: 5.4 COMPLEXITY
text: The time complexity of our Adaimp update for a batch of size B is O(D · d · B) which is similar, up to constants, to the complexity of updates for the AdaGrad algorithm. We do not notice any runtime differences between our algorithm applied in dimensions (d,D) and a CP decomposition of dimension D (see Section 7). The runtime for large enough D is dominated by the matrix product (O(D2 ·B)) required to compute the cross-entropy in Equation (1).
Algorithm 1 Step of PComplEx training (i, j, k)← sample from S gi, gj , gk ← gradients in (Pui, Puj , wk) ui, G̃i ← Adaimp (η, ui, gu, G̃i, P ) uj , G̃j ← Adaimp (η, uj , gj , G̃j , P ) wk, Gk ← AdaGrad(η, wk, gk, Gk)
Algorithm 2 Adaimp
Input: η, x(t), g(t), G̃(t−1), P g̃(t) ← P (P⊤g(t)) G̃(t) ← G̃(t−1) + g̃(t) ⊙ g̃(t)
x(t+1) ← x(t) − ηP⊤Diag(G̃(t))−1/2g̃(t) return x(t+1), G̃(t)
Figure 2: The two algorithms used in the training of PComplEx

heading: 6 PROJECTED COMPLEX
text: As the state-of-the-art variant of CP is ComplEx (Trouillon et al., 2016; Lacroix et al., 2018), we propose the following alternative to PCP base on ComplEx with Adaimp in practice. Given the ComplEx decomposition X = Re([U, V, U ]), a low-rank decomposition of the entity factor U as PŨ leads to the model PComplEx we use in the experiments of Section 7:
(PCOMPLEX)Xijk = ⟨Pui, vj , Puk⟩ = ⟨ui, P⊤Diag(vj)P, uk⟩, ui, wk ∈ Cd, vj ∈ CD, fixed P ∈ RD×d
PComplEx is similar to ComplEx but with interactions described by full matrices of rank D that share a same basis. We learn this decomposition with Algorithms 1 and 4.

heading: 7 EXPERIMENTS
text: In this section, we compare ComplEx optimized with AdaGrad and PComplEx optimized with Adaimp. We optimize the regularized empirical risk of Equation (2). Following Lacroix et al. (2018), we regularize ComplEx with the weighted nuclear-3 norm, which is equivalent to regularizing ∥ui∥33 + ∥uj∥33 + ∥wk∥33 for each training example (i, j, k). For PComplEx based models, we regularize ∥Pui∥33 + ∥vj∥33 + ∥Puk∥33 by analogy.
We conduct all experiments on a Quadro GP 100 GPU. The code for PComplEx and Adaimp is available in the supplementary materials, experiments on ComplEx use the code1 from (Lacroix et al., 2018). We include results from TuckER (Balažević et al., 2019b), DRT and SRT which are the two models considered in Wang et al. (2019), ConvE (Dettmers et al., 2018), HolEx (Xue et al., 2018), LFM (Jenatton et al., 2012) and MurP (Balažević et al., 2019a) without re-implementation on our parts. All the parameters for the experiments in this section are reported in Appendix 9.11.

heading: 7.1 DATASETS
text: WN18 (Bordes et al., 2013) is taken from the Wordnet database which contains words and relations between them. WN18RR (Dettmers et al., 2018) is a filtering of this dataset which removes train/test leakage. YAGO3-10 (Dettmers et al., 2018) is taken from the eponymous knowledge base. Finally, SVO (Jenatton et al., 2012) contains observations of Subject, Verb, Object triples. All statistics for these datasets can be found in Appendix 9.13. Experiments on the FB15K and FB15K-237 datasets are deferred to Appendix 9.10.

heading: 7.2 RESULTS
text: We report the filtered Mean Reciprocal Rank (Nickel et al., 2016b) on Figure 3. For SVO, we report the only figure available in previous work which is the filtered hits at 5% (Jenatton et al., 2012). These measures are detailed in Appendix 9.11. Only the grid-search parameters were given for LFM, so we were not able to obtain a precise number of parameters for the number they report.
1available at https://github.com/facebookresearch/kbc
On WN18, SVO and YAGO3-10 we observe sizable performance gains for low embedding sizes : up to 0.14 MRR points on WN18, 0.05 MRR points on YAGO and 0.03 H@5% points on SVO.
The TuckER (Balažević et al., 2019b) model performs similarly to PComplEx and ComplEx except on FB15K and WN18 where it underperforms (see Appendix 9.10). We expect this discrepancy to come from a less extensive grid-search rather than any intrinsic differences in the models that are both based on the Tucker decomposition. The consistency on all operating points of our method with ComplEx shows an advantage of our method, which enjoys the same learning rate robustness as AdaGrad, and does not require choosing a learning-rate decay, leading to easier experiments with only the regularization strength to tune. The MurP model (Balažević et al., 2019a) provides good performances for low embedding sizes on WN18RR, but underperforms on FB15K-237 (see Appendix 9.10). All other models fail to match the performances of ComplEx and PComplEx with equivalent number of parameters.
Variance of performances in PComplEx due to random choice of P is similar to the variance of ComplEx. We present experiments on the WN18 dataset for 5 different seeds in Appendix 9.9.

heading: 8 CONCLUSION
text: By observing that the core tensor of the Tucker decomposition can itself be decomposed, we obtain new models that are reminiscent of the canonical decomposition with low-rank factors. We provide experimental evidence that a naive application of AdaGrad on this decomposition fails, due to individual coordinates losing their meaning. We propose a new algorithm, Adaimp, which fixes this issue. Our model, when optimized with Adaimp, provides better performances than ComplEx in the low-rank regime, and matches its performances in the other regimes.

heading: 9 APPENDIX
text: In subsections 9.1 and 9.2 where we discuss the Adagrad algorithm, we do so in a general setting, where we study, for fixed P with orthonormal columns, the problem :
min θ f(Pθ).

heading: 9.1 PROJECTED ADAGRAD UPDATE
text: Let D denote Diag(G)1/2 for the classical version of Adagrad and G1/2 itself for the “full” version of the update. When the parameter θ is constrained to a set Θ, the update proposed in Eq.(1) in Duchi et al. (2011) the one obtained by solving
min θ∈Θ
(θ̃ − z)⊤D(θ̃ − z) with z = θ̃(t) − ηD−1g(t)
To enforce the constraint that θ̃ = Pθ, we can consider the Lagrangian
L(θ̃, θ;λ) = (θ̃ − z)⊤D(θ̃ − z)− λ⊤(θ̃ − Pθ)
whose stationary points satisfy D(θ̃ − z) = λ and P⊤λ = 0. So this entails P⊤D(θ̃ − θ̃(t)) = ηP⊤g(t) and finally using θ̃ = Pθ we obtain an update in θ as follows
θ(t+1) = θ(t) − η(P⊤DP⊤)−1Pg(t). Clearly, PDP⊤ is in general non-diagonal whether D is diagonal or not, and so this approach does not provide a computationally efficient update.
If D = G1/2, then since PG1/2P⊤ = (PGP⊤)1/2 the update is the same as the full Adagrad update (5) that we derive in the following section and replacing (PGP⊤)1/2 by its diagonal approximation recovers update (3).

heading: 9.2 THE TWO FULL ADAGRAD UPDATES AND THE QUALITY OF APPROXIMATIONS OF THE
text: DIAG VERSIONS
If we consider the full versions of the Adagrad updates then (letting again Π = PP⊤) its application to θ 7→ f(Pθ) and θ̃ 7→ f(Πθ̃) yield respectively
θ̃(t+1) = θ̃(t) − η P ( P⊤G(t)P )−1/2 P⊤g(t) and (5)
θ̃(t+1) = θ̃(t) − ηΠ ( ΠG(t) Π )−†/2 Π g(t), (6)
where M† notes the pseudo-inverse of a matrix M . As it turns out, the two updates are equivalent:
Indeed, first (ΠGΠ)1/2 = P (P⊤GP )1/2P⊤ because P⊤P = I implies that
P (P⊤GP )1/2P⊤P (P⊤GP )1/2P⊤ = ΠGΠ,
and the p.s.d. squareroot is unique. Second, taking the pseudo-inverse of this identity, we have (ΠGΠ)†/2 = ( P (P⊤GP )1/2P⊤ )† = P (P⊤GP )−1/2P⊤. (7)
because, if H = (P⊤GP )1/2 is an invertible matrix, then (PHP⊤)† = PH−1P⊤ given that PHP⊤ PH−1P⊤ = PP⊤. Finally multiplying both sides of Eq. (7) by Π shows that
Π(ΠG(t)Π)†/2Π = P (P⊤G(t)P )−1/2P⊤.
This shows that although Adagrad is not invariant for any invertible P , it is invariant for any P such that P⊤P = I . Eq.(5) seems in general simpler than (6), but note that if D = d, then Π is the identity and (6) shows that both full updates are in that case actually equivalent to the full update of Adagrad applied to plain CP.
Finally, the equivalence of the full updates discussed above strongly suggests that if our proposed update performs better than the naive application of Adagrad to PCP, it is because PDiag(ΠG(t)Π)−1/2P⊤ is a better approximation of (G(t))−1/2 than Diag(PG(t)P⊤)−1/2 while not being much more computationally expensive.

heading: 9.3 CP-TUCKER
text: We have X = [[C;U, V,W ]] and C = [[P1, P2, P3]]. For all i, j, k, we have:
Xi,j,k = d∑ r1,r2,r3 Cr1,r2,r3Ui,r1Vi,r2Wk,r3
= d∑ r1,r2,r3
( D∑
s=0
[P1]r1,s[P2]r2,s[P3]r3,s ) Ui,r1Vi,r2Wk,r3
= D∑ s=0 ( d∑ r1 Ui,r1 [P1]r1,s )( d∑ r2 Vj,r2 [P2]r2,s )( d∑ r3 Wk,r3 [P3]r3,s ) = ⟨P1ui, P2vj , P3wk⟩

heading: 9.4 HOLEX AND LATENT FACTOR MODEL

heading: 9.4.1 HOLEX
text: The HolEx model (Xue et al., 2018) writes Xi,j,k = ∑R
r=1⟨(cr⊙ui) ⋆uj , wrk⟩, where ⋆ denotes the circular correlation2. Exploiting the equivalence between HolE and ComplEx (Hayashi & Shimbo, 2017), denoting byF the discrete Fourier transform (with values in C), we can write for embeddings of size d:
R∑ r=1 ⟨(cr ⊙ ui) ⋆ uj , wrk⟩ = 1 d Re
( R∑
r=1
⟨F(cr ⊙ ui),F(uj),F(wrk)⟩
)
= 1
d Re
( R∑
r=1
⟨F(cr) ⋆ F(ui),F(uj),F(wrk)⟩
)
For all vectors, we write ûi = F(ui) ∈ Cd. We can re-write the circular correlation F(cr) ⋆ F(ui) as Crûi where Cr ∈ Cd×d is the circulant matrix associated with cr ∈ Rd. We have :
R∑ r=1 ⟨(cr ⊙ ui) ⋆ vj , wrk⟩ = 1 d Re
( R∑
r=1
⟨Crûi, ûj , ŵrk⟩
) .
Finally, with C1 = [C1, . . . , CR] ∈ RRd×d the vertical stacking of all Cr, C2 = [Id, . . . , Id] and ŵk = [ŵ 1 k, . . . , ŵ R k ]:
R∑ r=1 ⟨(cr ⊙ ui) ⋆ uj , wrk⟩ = 1 d Re ( ⟨C1ûi, C2ûj , ŵk⟩ ) HolEx with embeddings of size d is close to the CPT with D = Rd, allowing for two different complex matrices to act on left and right hand side embeddings.

heading: 9.4.2 LATENT FACTOR MODEL
text: The latent factor model defines the score for a triple (i, j, k) as:
Xi,j,k = ⟨(si + z), Rj(ok + z′)⟩, with Rj = D∑
r=1
αjrurv ⊤ r .
Removing the bias terms z and z′ and gathering ur and vr into matrices P1 and P3 leads to the model CPT. In the PCP model, we fix P1 and P3 instead of learning them. We do not use a sparsity inducing penalty on α but rather a variational form of the nuclear norm on the whole tensor.
2[a ⋆ b]k = ∑d−1
i=0 aib(i+k) mod d.

heading: 9.5 TUCKER2 WITH CP-TUCKER
text: Let for all 1 ≤ r ≤ d, M (r) be a matrix of zeros except its r− th column which is all one. Let P1 be the vertical concatenation of all (M (r))r=1..d and P2 the vertical concatenation of d identity matrix in Rd. Remember that for all k, wk is an element of Rd 2
. For all 0 ≤ r < d, let wrk be the restriction of wk to its [rd, (r + 1)d] coordinates.
Then P1 and P2 are elements of Rd 2×d and we have for all i, j, k:
⟨P1ui, P2vj , wk⟩ = d∑
r1=1
⟨M (r1)ui, Idvj , wr1k ⟩
= d∑ r1=1 ui,r1⟨vj , w r1 k ⟩ by definition of M (r1)
= d∑ r1=1 ui,r1
( d∑
r2=1
vj,r2wk,r1r2 ) by definition of wr1k
= u⊤i Mat(wk)vj

heading: 9.6 COMPLETE PCOMPLEX ALGORITHM
text: Let U ∈ CN×d be the the entity embedding matrix and V ∈ C2L×D be the predicate embedding matrix. Let P : RN×d 7→ RN×D such that P(Ui) = PUi. Let GUt and GVt be the stochastic gradients with respect to U and V at time t.
Algorithm 3 PComplEx optimized with Adaimp
Input: learning rate η, (random) matrix P with orthogonal columns, ϵ CU , CV ← 0 while U, V not converged do G̃Ut ← P(GUt ) CU ← CU + G̃Ut ⊙ G̃Ut CV ← CV +GVt ⊙GVt U ← U − η · P⊤(G̃Ut /( √ CU + ϵ))
V ← V − η ·GVt /( √ CV + ϵ))
end while return U

heading: 9.7 ADAM - IMPLICIT
text: For WN18RR, we used the ideas presented in this paper, but adapted to the Adam (Kingma & Ba, 2014) optimization algorithm. Similar to Adaimp, first and second moment estimates are accumulated in RD and we project back to Rd only for the parameter update. For simplicity, we present here the dense version of the algorithm applied to the entity embeddings U ∈ RN×d. Let P : RN×d 7→ RN×D such that P(Ui) = PUi. Let Gt be the stochastic gradient with respect to U at time t.
9.8 ADArow
For the same control experiment as in Section 5, we observe the performances of Adarow which is rotation invariant by design.
Algorithm 4 Adamimp applied to U Input: η, β1, β2,P, ϵ m0, v0, t← 0 while U not converged do t← t+ 1 G̃t ← P(Gt) mt ← β1 ·mt−1 + (1− β1) · G̃t vt ← β2 · vt−1 + (1− β2) · G̃t ⊙ G̃t m̂t ← mt/(1− βt1) v̂t ← vt/(1− βt2) U ← U − η · P⊤(m̂t/( √ v̂t + ϵ))
end while return U

heading: 9.9 VARIANCE OF PCOMPLEX
text: We run 5 grid search and plot the 5 associated convex-hulls on the WN18 dataset optimized with Adaimp. Note that despite the added randomness of the choice of P , there is not more variance in PComplEx than in ComplEx.
0 20 40 60 80 100
params per entities
0.2
0.4
0.6
0.8
M R R
Variance of PComplEx on WN18
ComplEx
PComplEx

heading: 9.10 FB15K DATASETS
text: We use two subsets of the Freebase knowledge base : FB15K (Bordes et al., 2013) and FB15K237 (Toutanova & Chen, 2015). FB15K-237 is a harder version of FB15K, where some triples have been removed to avoid leakage between the train and test set. There is no difference of performances between PComplEx and ComplEx on these datasets.
10 1
10 2
10 3
Params per entities
0.4
0.6
0.8
M R R
FB15K
ComplEx
PComplEx
TuckER
HolEx
10 1
10 2
10 3
10 4
Params per entities
0.275
0.300
0.325
0.350
M R R
FB15K-237
TuckER
D/SRT
MurP
ConvE

heading: 9.11 EXPERIMENTAL DETAILS
text: Metrics Let rank(X̂i,j,:; k) be the rank of X̂i,j,k in the sorted list of values of X̂i,j,:. We report the MRR for most datasets :
MRR(X) = 1 |S| ∑
(i,j,k)∈S
1
rank(Xi,j,:; k) .
For SVO, the task is slightly different as the ranking happens on the predicate mode. The metric reported in Jenatton et al. (2012) is Hits@5% defined as :
H@5%(X) = 1 |S| ∑
(i,j,k)∈S
1(rank(Xi,:,k; j) ≤ 227).
The metrics we report are filtered by ignoring other true positives when computing the ranks, as done in Bordes et al. (2013).
Number of parameters per entities We count the number of floats in the model, and divide by the number of entities in the dataset. For different methods, the number of parameters are :
• ComplEx: 2 · d · (N + 2L) • PComplEx: 2 · d ·N + 2 ·D · 2L+ d ·D • TuckEr: N · de + L · dp + d2e · dp • MurP: N · (d+ 1) + 2 · L · d • ConvE: taken from Dettmers et al. (2018) • D/SRT: taken from Wang et al. (2019) • HolEx: d · (N + L)
Grid Search For SVO:
• For Complex vary d in [5, 25, 50, 100, 500, 1000, 2000]. For PComplEx, we vary d in [5, 25, 50, 100, 500]. • The strength of the regularizer, ν varies in [5e− 4, 1e− 3, 5e− 3, 1e− 2, 5e− 2, 1e− 1]. • Finally, for PComplEx, we vary the dimension D in [5, 25, 50, 100, 500, 1000, 2000, 4000, 8000].
For all other datasets, we run 500 epochs :
• For FB15K and FB15K-237, we vary d in [5, 25, 50, 100, 500, 1000, 2000]. For YAGO3-10, WN18 and WN18RR, we add ranks 8 and 16 to that list. • The strength of the regularizer, ν varies in [5e− 4, 1e− 3, 5e− 3, 1e− 2, 5e− 2, 1e− 1]. • Finally, for PComplEx, we vary the dimension D in [5, 25, 50, 100, 500, 1000, 2000].
For TuckEr, we use the hyperparameters described in Balažević et al. (2019b) for each dataset. We apply a multiplicative factor γ to the dropout rates and run a grid-search over this parameter.
• For WN18RR, we vary de in [12, 15, 18, 20, 40, 100, 200]. For FB15K-237, we vary de in [10, 20, 30, 40, 50, 80, 100, 200]. For WN18, de varies in [5, 8, 16, 25, 50, 100]. For FB15K de varies in . • For WN18RR, WN18 we search γ over [0, 0.5, 1]. For FB15K-237 where this lead to non
increasing performances as a function of number of parameters, we refined this grid to [0, 0.1, 0.2, ..., 0.7] for de ≤ 80. For FB15K, we use a grid of [0, 0.2, 0.4, ...1].
For MurP, we use the hyperparameters described in Balažević et al. (2019a) for each dataset and vary d in [10, 12, 15, 18, 20, 40, 100, 200] on WN18RR and in [10, 20, 40, 100, 200] for FB15K-237.
Other details Batch-size is fixed to 1000, the learning-rate is 1e − 1 both for Adagrad and for Adaimp, we run 100 epochs to ensure convergence. The best model based on validation MRR is selected and we report the corresponding test MRR.

heading: 9.12 RUNNING TIMES
text: We give here the running times of each method per epoch, on the WN18RR dataset for comparable dimensionalities and run on a P100 GPU. We use the original implementation for MurP (Balažević et al., 2019a) and TuckEr (Balažević et al., 2019b). For ComplEx, we use the implementation from Lacroix et al. (2018).
• ComplEx (d = 200) : 4s/epoch. • PComplEx (d = D = 200, Adaimp) : 5s/epoch. • MurP (d = 200): 38s/epoch. • TuckEr (de = 200, dr = 200): 24s/epoch
Note that these running times are implementation dependent. In the figure below, we give the learning-curves of MurP, TuckEr, ComplEx and PComplEx for one operating point on the WN18RR dataset. The convergence speed of these methods is given in the figure below at an operating point on WN18RR where all methods are close in final performances.
100 200 300 400 500
Epoch
0.35
0.40
0.45
0.50
M R R
Convergence on WN18RR
MurP
TuckEr
ComplEx
PComplEx

heading: 9.13 DATASET STATISTICS

Review summaries:
In this paper, a tensor decomposition method is studied for link prediction problems. The model is based on Tucker decomposition but the core tensor is decomposed as CP decomposition so that it can be seen as an interpolation between Tucker and CP. The performance is evaluated with several NLP data sets (e.g., subject-verb-object triplets). 

The paper introduces a novel tensor decomposition that is reminiscent of canonical decomposition (CP) with low-rank factors, based on the observation that the core tensor in Tucker decomposition can be decomposed, resulting in a model interpolating between CP and Tucker. The authors argue that a straight application of AdaGrad on this decomposition is inadequate, and propose Ada^{imp} algorithm that enforces rotation invariance of the gradient update. The new decomposition is applied to ComplEx model (called PComplEx) that demonstrates better performance than the baseline.