Title: PROJECTED CANONICAL DECOMPOSITION FOR KNOWLEDGE BASE COMPLETION

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.

1. INTRODUCTION:

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žević 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žević 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žević 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žević 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.

2. LINK PREDICTION IN KNOWLEDGE BASES:

Notation Tensors and matrices are denoted by uppercase letters. For a matrix U , u i 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 ⊙.

2.1. LEARNING SETUP:

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) = -X i,j,k + log ( ∑ k ′ ̸ =k exp(X i,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.

2.2. RELATED WORK:



2.2.1. CANONICAL DECOMPOSITION AND ITS VARIANTS:

The canonical polyadic (CP) decomposition of a tensor X is defined entrywise by

∀i, j, k, X i,j,k = ⟨u i , v j , w k ⟩ := d ∑ r=1 u ir v jr w kr .

The smallest value of d for which this decomposition exists is the rank of X. Each element X i,j,k is thus represented as a multi-linear product of the 3 embeddings in R d 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., w k = ūk . Prior to ComplEx, Dismult was proposed (Yang et al., 2014) as a variant of CP with w k = u k . 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.

2.2.2. TUCKER DECOMPOSITION AND ITS VARIANTS:

Given a tensor X of size N ×L×N , the Tucker decomposition of X is defined entrywise by

∀i, j, k, X i,j,k = ⟨u i ⊗ v j ⊗ w k , C⟩ := d1 ∑ r1=1 d2 ∑ r2=1 d3 ∑ r3=1 C r1,r2,r3 u ir1 v jr2 w kr3 .

The triple (d 1 , d 2 , d 3 ) 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 u j , v k , w l and C the three dimensional d 1 × d 2 × d 3 core tensor. Note that the CP decomposition is a Tucker decomposition in which

d 1 = d 2 = d 3 = 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žević 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.

2.2.3. OTHER APPROACHES:

(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žević 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.

3. INTERPOLATING BETWEEN CP AND TUCKER:

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):

X ijk = ⟨u i ⊗ v j ⊗ w k , C⟩ = ⟨P 1 u i , P 2 v j , P 3 w k ⟩ or equivalently X = [[U P ⊤ 1 , V P ⊤ 2 , W P ⊤ 3 ]],

where P 1 , P 2 , P 3 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 P 2 = I. Learning matrices U, V, W, P 1 , P 3 of this decomposition simultaneously leads to the following model, which we call CP-Tucker (CPT):

(CPT) X ijk = ⟨P 1 u i , v j , P 3 w k ⟩, u i , w k ∈ R d , v j ∈ R D , P i ∈ R D×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 P 1 and P 3 are full rank), whereas for D = d 2 we recover a full Tucker model because the matrices P 1 and P 3 can be chosen such that

⟨P 1 u i , v j , P 3 w k ⟩ = u i Mat(v j )v T j ,

where Mat is the operator that maps a d 2 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 P i . 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 P 1 and P 3 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 P 1 , P 3 project the embeddings of dimension d into a higher dimension D: 

(PCP) X ijk = ⟨P 1 u i , v j , P 3 w k ⟩, u i , w k ∈ R d , v j ∈ R D , fixed P 1 , P 3 ∈ R D×d

4. MOTIVATION: OPTIMIZATION ISSUES WITH CPT AND PCP:

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.

4.1 CONTROL EXPERIMENT: UNITARY P 1 AND P 3 IN PCP When D = d in PCP, the model becomes equivalent to CP. Indeed, the matrices P 1 and P 3 are unitary (P

1 P ⊤ 1 = P 3 P ⊤ 3 = I) and so [[(U P 1 )P ⊤ 1 , V, (W P 3 )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 P 1 and P 2 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 P 1 and P 2 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 P 1 and P 3 (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: ADA imp:

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.

5.1. TWO EQUIVALENT PARAMETRIZATIONS OF PCP:

The following decomposition is equivalent to PCP, but its embeddings are expressed in R D :

(PCP FULL ) X i,j,k = ⟨P 1 P ⊤ 1 u i , v j , P 3 P ⊤ 3 w k ⟩, with u i , v j , w k ∈ R D , fixed P 1 , P 2 ∈ R D×d .

Note that with u i = P ⊤ 1 u i and w k = P ⊤ 3 w k , we go from PCP full to PCP. The practical differences between PCP and PCP full are that PCP full learn embeddings in the high-dimensional space, maintaining the low-rank structure of the overall entity embeddings through the orthogonal projections P 1 P T 1 and P 3 P T 3 . The practical interest of PCP full 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, P 1 and P 3 disappear in PCP full , so that optimizing a PCP full 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: Ada imp The approach we propose is to effectively optimize PCP full with Adagrad. However, when d < D, which is the interesting case for PCP, we notice that we do not need to maintain embeddings in R D . Our approach, called Ada imp , computes the gradients and Adagrad coefficients with respect to u i , w k ∈ R D , but the full dimension factor matrices U and W are never explicitly stored in memory. Rather, we store

u i = P ⊤ 1 u i and w k = P ⊤ 3 w k ∈ R d ,

which is all that is required for any model computation in PCP full since P 1 and P 3 are fixed. Overall, the effective model parameters are exactly the same as in PCP, and we call this approach PCP +Ada imp .

An Ada imp update is described in Algorithm 4. While PCP +Adagrad and PCP +Ada imp work with the same number of model parameters, the fundamental difference is the computation of Adagrad coefficients. Since Ada imp effectively applies Adagrad to PCP full , we need to maintain the Adagrad coefficients in R D even when d < D: the overall update is first computed in R D and projected back to R d after the application of Adagrad rescaling. In constrat, in vanilla PCP +Adagrad, the gradient is projected to R d before Adagrad rescaling.

5.2. IMPLICIT OPTIMIZATION OF PCP FULL : ADA imp:

In this section, we discuss more formally the Ada imp updates, and how they compare to PCP +Adagrad. In the following, Ũ , W are in R N ×d , whereas U, V and W are in R N ×D . Using d ≤ D, P 1 , P 3 ∈ R D×d , and we use the notation Π 1 = P 1 P ⊤ 1 and Π 3 = P 3 P ⊤ 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 Ũ ,V, W L( Ũ P ⊤ 1 , V, W P ⊤ 3 ) (PCP FULL ) argmin U,V,W L(U Π ⊤ 1 , V, W Π ⊤ 3 )

We focus on a step at time (t) on vectors ũi and u i . We assume that at this time t, the tensor iterates are the same, that is

Ũ = U P ⊤ 1 (resp. W = W P ⊤ 1 ) so that [[ Ũ 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 = ϵI d + ∑ t τ =1 g (τ ) i g (τ ) i ⊤ .

The updates for (PCP) are:

ũt+1 i = ũt i -ηDiag(P ⊤ 1 G (t) i P 1 ) -1/2 P ⊤ 1 g (t)

i .

(3)

Note that due to the presence of P 1 inside the Diag operator, the update (3) is not rotation invariant. Moreover, for random P 1 , the matrix P ⊤ 1 G

(t)

i P 1 will be far from diagonal with high probability, making its diagonal meaningless. This is visualized in Figure 1b .

Similar updates can be derived for (PCP full ):

u t+1 i = 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 ∈ R N ×D , we apply these updates implicitly with Ada imp , by storing u

(t) i = P ⊤ 1 u (t) i in R d .

Let us compare the updates :

(ADAGRAD) ũt+1 i = ũt i -η × Adagrad in R d Diag(P ⊤ 1 G (t) i P 1 ) -1/2 × projection to R d P ⊤ 1 g (t) i (ADA imp ) u t+1 i = u t i -η × P ⊤ 1 projection to R d × Diag(Π ⊤ 1 G (t) i Π 1 ) -1/2 Adagrad in R D × Π ⊤ 1 g (t) i projection to Im(Π)∈R D

Going back to our control experiment, we note on Figure 1a that PCP +Ada imp matches the performances of CP+Adagrad for all D, indicating that we fixed this optimization issue.

5.3. ALTERNATIVES TO ADA imp:

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 Ada row ) 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 Ada imp are indeed better is further backed by the experiments in the next section .

5.4. COMPLEXITY:

The time complexity of our Ada imp 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(D 2 • B)) required to compute the cross-entropy in Equation (1).

Algorithm 1 Step of PComplEx training:

(i, j, k) ← sample from S g i , g j , g k ← gradients in (P u i , P u j , w k ) u i , Gi ← Ada imp (η, u i , g u , Gi , P ) u j , Gj ← Ada imp (η, u j , g j , Gj , P ) w k , G k ← AdaGrad(η, w k , g k , G k )

Algorithm 2 Ada imp 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/2 g(t) return x (t+1) , G(t)

6. PROJECTED COMPLEX:

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 Ada imp in practice. Given the ComplEx decomposition X = Re([U, V, U ]), a low-rank decomposition of the entity factor U as P Ũ leads to the model PComplEx we use in the experiments of Section 7:

(PCOMPLEX) X ijk = ⟨P u i , v j , P u k ⟩ = ⟨u i , P ⊤ Diag(v j )P, u k ⟩, u i , w k ∈ C d , v j ∈ C D , fixed P ∈ R D×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.

7. EXPERIMENTS:

In this section, we compare ComplEx optimized with AdaGrad and PComplEx optimized with Ada imp . 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 ∥u i ∥ 3 3 + ∥u j ∥ 3 3 + ∥w k ∥ 3 3 for each training example (i, j, k). For PComplEx based models, we regularize ∥P u i ∥ 3 3 + ∥v j ∥ 3 3 + ∥P u k ∥ 3 3 by analogy. We conduct all experiments on a Quadro GP 100 GPU. The code for PComplEx and Ada imp is available in the supplementary materials, experiments on ComplEx use the code1 from (Lacroix et al., 2018) . We include results from TuckER (Balažević 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žević et al., 2019a) without re-implementation on our parts. All the parameters for the experiments in this section are reported in Appendix 9.11.

7.1 DATASETS 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.

7.2. RESULTS:

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. For some datasets (WN18, WN18RR), Adam (Kingma & Ba, 2014) is beneficial, in which case we use the implicit adaptation of Adam detailed in Appendix 9.7.

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žević 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žević 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.

8. CONCLUSION:

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, Ada imp , which fixes this issue. Our model, when optimized with Ada imp , provides better performances than ComplEx in the low-rank regime, and matches its performances in the other regimes.

9. APPENDIX:

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 θ).

9.1. PROJECTED ADAGRAD UPDATE:

Let D denote Diag(G) 1/2 for the classical version of Adagrad and G 1/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 -1 g (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 ⊤ ) -1 P g (t) .

Clearly, P DP ⊤ is in general non-diagonal whether D is diagonal or not, and so this approach does not provide a computationally efficient update.

If D = G 1/2 , then since P G 1/2 P ⊤ = (P GP ⊤ ) 1/2 the update is the same as the full Adagrad update (5) that we derive in the following section and replacing (P GP ⊤ ) 1/2 by its diagonal approximation recovers update (3).

9.2. THE TWO FULL ADAGRAD UPDATES AND THE QUALITY OF APPROXIMATIONS OF THE DIAG VERSIONS:

If we consider the full versions of the Adagrad updates then (letting again Π = P P ⊤ ) its application to θ → f (P θ) and θ → f (Π θ) yield respectively θ(t+1) = θ(t) -η P ( P ⊤ G (t) P ) -1/2 P ⊤ g (t) and ( 5)

EQUATION

)

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/2 P ⊤ because P ⊤ P = I implies that P (P ⊤ GP ) 1/2 P ⊤ P (P ⊤ GP ) 1/2 P ⊤ = ΠGΠ, and the p.s.d. squareroot is unique. Second, taking the pseudo-inverse of this identity, we have

EQUATION

)

because, if H = (P ⊤ GP ) 1/2 is an invertible matrix, then (P HP ⊤ ) † = P H -1 P ⊤ given that P HP ⊤ P H -1 P ⊤ = P P ⊤ . Finally multiplying both sides of Eq. ( 7) by Π shows that Π(ΠG (t) Π) †/2 Π = P (P ⊤ G (t) P ) -1/2 P ⊤ .

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 P Diag(ΠG (t) Π) -1/2 P ⊤ is a better approximation of (G (t) ) -1/2 than Diag(P G (t) P ⊤ ) -1/2 while not being much more computationally expensive.

9.3. CP-TUCKER:

We have X = [[C; U, V, W ]] and C = [[P 1 , P 2 , P 3 ]]. For all i, j, k, we have:

X i,j,k = d ∑ r1,r2,r3 C r1,r2,r3 U i,r1 V i,r2 W k,r3 = d ∑ r1,r2,r3 ( D ∑ s=0 [P 1 ] r1,s [P 2 ] r2,s [P 3 ] r3,s ) U i,r1 V i,r2 W k,r3 = D ∑ s=0 ( d ∑ r1 U i,r1 [P 1 ] r1,s ) ( d ∑ r2 V j,r2 [P 2 ] r2,s ) ( d ∑ r3 W k,r3 [P 3 ] r3,s ) = ⟨P 1 u i , P 2 v j , P 3 w k ⟩

9.4 HOLEX AND LATENT FACTOR MODEL

9.4.1. HOLEX:

The HolEx model (Xue et al., 2018) writes

X i,j,k = ∑ R r=1 ⟨(c r ⊙ u i ) ⋆ u j , w r k ⟩,

where ⋆ denotes the circular correlation2 . Exploiting the equivalence between HolE and ComplEx (Hayashi & Shimbo, 2017) , denoting by F the discrete Fourier transform (with values in C), we can write for embeddings of size d:

R ∑ r=1 ⟨(c r ⊙ u i ) ⋆ u j , w r k ⟩ = 1 d Re ( R ∑ r=1 ⟨F(c r ⊙ u i ), F(u j ), F(w r k )⟩ ) = 1 d Re ( R ∑ r=1 ⟨F(c r ) ⋆ F(u i ), F(u j ), F(w r k )⟩ )

For all vectors, we write ûi = F(u i ) ∈ C d . We can re-write the circular correlation F(c r ) ⋆ F(u i ) as C r ûi where C r ∈ C d×d is the circulant matrix associated with c r ∈ R d . We have :

R ∑ r=1 ⟨(c r ⊙ u i ) ⋆ v j , w r k ⟩ = 1 d Re ( R ∑ r=1 ⟨C r ûi , ûj , ŵr k ⟩ ) . Finally, with C 1 = [C 1 , . . . , C R ] ∈ R Rd×d the vertical stacking of all C r , C 2 = [I d , . . . , I d ] and ŵk = [ ŵ1 k , . . . , ŵR k ]: R ∑ r=1 ⟨(c r ⊙ u i ) ⋆ u j , w r k ⟩ = 1 d Re ( ⟨C 1 ûi , C 2 ûj , ŵ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.

9.4.2. LATENT FACTOR MODEL:

The latent factor model defines the score for a triple (i, j, k) as:

X i,j,k = ⟨(s i + z), R j (o k + z ′ )⟩, with R j = D ∑ r=1 α j r u r v ⊤ r .

Removing the bias terms z and z ′ and gathering u r and v r into matrices P 1 and P 3 leads to the model CPT. In the PCP model, we fix P 1 and P 3 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.

9.5. TUCKER2 WITH CP-TUCKER:

Let for all 1 ≤ r ≤ d, M (r) be a matrix of zeros except its r -th column which is all one. Let P 1 be the vertical concatenation of all (M (r) ) r=1..d and P 2 the vertical concatenation of d identity matrix in R d . Remember that for all k, w k is an element of R d 2 . For all 0 ≤ r < d, let w r k be the restriction of w k to its [rd, (r + 1)d] coordinates.

Then P 1 and P 2 are elements of R d 2 ×d and we have for all i, j, k:

⟨P 1 u i , P 2 v j , w k ⟩ = d ∑ r1=1 ⟨M (r1) u i , I d v j , w r1 k ⟩ = d ∑ r1=1 u i,r1 ⟨v j , w r1 k ⟩ by definition of M (r1) = d ∑ r1=1 u i,r1 ( d ∑ r2=1 v j,r2 w k,r1r2 ) by definition of w r1 k = u ⊤ i Mat(w k )v j 9.6 COMPLETE PCOMPLEX ALGORITHM

Let U ∈ C N ×d be the the entity embedding matrix and V ∈ C 2L×D be the predicate embedding matrix. Let P : R N ×d → R N ×D such that P(U i ) = P U i . Let G U t and G V t be the stochastic gradients with respect to U and V at time t.

Algorithm 3 PComplEx optimized with Ada imp:

Input: learning rate η, (random) matrix P with orthogonal columns, ϵ

C U , C V ← 0 while U, V not converged do GU t ← P(G U t ) C U ← C U + GU t ⊙ GU t C V ← C V + G V t ⊙ G V t U ← U -η • P ⊤ ( GU t /( √ C U + ϵ)) V ← V -η • G V t /( √ C V + ϵ)) end while return U 9.7 ADAM -IMPLICIT

For WN18RR, we used the ideas presented in this paper, but adapted to the Adam (Kingma & Ba, 2014) optimization algorithm. Similar to Ada imp , first and second moment estimates are accumulated in R D and we project back to R d only for the parameter update. For simplicity, we present here the dense version of the algorithm applied to the entity embeddings U ∈ R N ×d . Let P : R N ×d → R N ×D such that P(U i ) = P U i . Let G t be the stochastic gradient with respect to U at time t.

9.8. ADA row:

For the same control experiment as in Section 5, we observe the performances of Ada row which is rotation invariant by design.

Algorithm 4 Adam imp applied to U 

Input: η, β 1 , β 2 , P, ϵ m 0 , v 0 , t ← 0 while U not converged do t ← t + 1 Gt ← P(G t ) m t ← β 1 • m t-1 + (1 -β 1 ) • Gt v t ← β 2 • v t-1 + (1 -β 2 ) • Gt ⊙ Gt mt ← m t /(1 -β t 1 ) vt ← v t /(1 -β t 2 ) U ← U -η • P ⊤ ( mt /( √ vt + ϵ)) end while return U

9.9. VARIANCE OF PCOMPLEX:

We run 5 grid search and plot the 5 associated convex-hulls on the WN18 dataset optimized with Ada imp . Note that despite the added randomness of the choice of P , there is not more variance in PComplEx than in ComplEx. We use two subsets of the Freebase knowledge base : FB15K (Bordes et al., 2013) and FB15K-237 (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. 

9.11. EXPERIMENTAL DETAILS:

Metrics Let rank( Xi,j,: ; k) be the rank of Xi,j,k in the sorted list of values of Xi,j,: . We report the MRR for most datasets :

M RR(X) = 1 |S| ∑ (i,j,k)∈S

1 rank(X i,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

(rank(X i,:,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 : • 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]. [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žević 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 d e in [12, 15, 18, 20, 40, 100, 200] . For FB15K-237, we vary d e in [10, 20, 30, 40, 50, 80, 100, 200] . For WN18, d e varies in [5, 8, 16, 25, 50, 100] . For FB15K d e 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 d e ≤ 80. For FB15K, we use a grid of [0, 0.2, 0.4, ...1].

For MurP, we use the hyperparameters described in Balažević 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 Ada imp , we run 100 epochs to ensure convergence. The best model based on validation MRR is selected and we report the corresponding test MRR.

9.12. RUNNING TIMES:

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žević et al., 2019a) and TuckEr (Balažević et al., 2019b) . For ComplEx, we use the implementation from Lacroix et al. (2018) . 

:

available at https://github.com/facebookresearch/kbc

[a ⋆ b] k = ∑ d-1 i=0 aib (i+k) mod d .

Figures and Tables:

imp ) (a) Performances on FB15K-237 in the control experiments for D = d, i.e., when PCP is a reparameterization of CP. We observe that PCP and CPT with vanilla Adagrad, which are variants of Tucker, underperform compared to CP. As expected in this case D = d, our modification of the Adagrad update leads to the same performances for PCP as for CP. Adagrad coefficients for subject/object embedding matrices in the control experiments (D = d) for CP and PCP, averaged by columns (i.e., embedding dimension) and sorted by values. Adagrad coefficients decay exponentially for CP, but the values are similar across most dimensions in PCP: the fixed unitary transform in PCP removes the benefit of Adagrad.

Figure 1: Results of the control experiment of Section 4.

When D = d the matrices P i 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 d 2 .

Figure 2: The two algorithms used in the training of PComplEx

Figure 3: MRR as a function of #floats / entities (see Appendix 9.11) on four knowledge bases. We plot the convex envelope of various operating points we tested, varying D for several values of d.For some datasets (WN18, WN18RR), Adam (Kingma & Ba, 2014) is beneficial, in which case we use the implicit adaptation of Adam detailed in Appendix 9.7.

Figure 4: Ada row vs AdaGrad on FB15K-237.

ComplEx: 2 • d • (N + 2L) • PComplEx: 2 • d • N + 2 • D • 2L + d • D • TuckEr: N • d e + L • d p + d 2 e • d p • 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:

ComplEx (d = 200) : 4s/epoch. • PComplEx (d = D = 200, Ada imp ) : 5s/epoch. • MurP (d = 200): 38s/epoch. • TuckEr (d e = 200, d r = 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.

Peer review (iclr2020):
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). Although the entire idea is interesting, the current form of the paper is not sufficient for acceptance. The main reasons are (A) the proposed model is not completely novel and (B) the empirical results are not significant. (A) The idea of combining CP and Tucker is not new. For example, Tomioka et al. (2010; Section 3.4) considered the Tucker-CP patterns (CP decomposition of the Tucker core). Although they used the Tucker-CP model to improve the interpretability rather than link prediction, the paper needs to make some attribution to the prior work. (B) By looking Figure 3, the proposed method, PComplEx, is not significantly better than the existing methods such as ComplEx. Except SVO data, PComplEx and ComplEx share almost the same performance curve. Also, other existing methods such as TuckER and MurP are evaluated only in a few points while (P)ComplEx is evaluated in many points. I feel this is unfair. Tomioka, R., Hayashi, K., & Kashima, H. (2010). Estimation of low-rank tensors via convex optimization. arXiv preprint arXiv:1010.0789.