Knowledge Graph Completion via Complex Tensor Factorization

Introduction

Web-scale knowledge graph provide a structured representation of world knowledge, with projects such as the Google Knowledge Vault Dong et al. (2014). They enable a wide range of applications including recommender systems, question answering and automated personal agents. The incompleteness of these knowledge graphs—also called knowledge bases—has stimulated research into predicting missing entries, a task known as link prediction or knowledge graph completion. The need for high quality predictions required by link prediction applications made it progressively become the main problem in statistical relational learning Getoor and Taskar (2007), a research field interested in relational data representation and modeling.

Knowledge graphs were born with the advent of the Semantic Web, pushed by the World Wide Web Consortium (W3C) recommendations. Namely, the Resource Description Framework (RDF) standard, that underlies knowledge graphs’ data representation, provides for the first time a common framework across all connected information systems to share their data under the same paradigm. Being more expressive than classical relational databases, all existing relational data can be translated into RDF knowledge graphs Sahoo et al. (2009).

Knowledge graphs express data as a directed graph with labeled edges (relations) between nodes (entities). Natural redundancies between the recorded relations often make it possible to fill in the missing entries of a knowledge graph. As an example, the relation CountryOfBirth could not be recorded for all entities, but it can be inferred if the relation CityOfBirth is known. The goal of link prediction is the automatic discovery of such regularities. However, many relations are non-deterministic: the combination of the two facts IsBornIn(John,Athens) and IsLocatedIn(Athens,Greece) does not always imply the fact HasNationality(John,Greece). Hence, it is natural to handle inference probabilistically, and jointly with other facts involving these relations and entities. To this end, an increasingly popular method is to state the knowledge graph completion task as a 3D binary tensor completion problem, where each tensor slice is the adjacency matrix of one relation in the knowledge graph, and compute a decomposition of this partially-observed tensor from which its missing entries can be completed.

Factorization models with low-rank embeddings were popularized by the Netflix challenge Koren et al. (2009). A partially-observed matrix or tensor is decomposed into a product of embedding matrices with much smaller dimensions, resulting in fixed-dimensional vector representations for each entity and relation in the graph, that allow completion of the missing entries. For a given fact r(s,o) in which the subject entity $s$ is linked to the object entity $o$ through the relation $r$, a score for the fact can be recovered as a multilinear product between the embedding vectors of $s$, $r$ and $o$, or through more sophisticated composition functions Nickel et al. (2016a).

Binary relations in knowledge graphs exhibit various types of patterns: hierarchies and compositions like FatherOf, OlderThan or IsPartOf, with strict/non-strict orders or preorders, and equivalence relations like IsSimilarTo. These characteristics maps to different combinations of the following properties: reflexivity/irreflexivity, symmetry/antisymmetry and transitivity. As described in Bordes et al. (2013a), a relational model should (i) be able to learn all combinations of such properties, and (ii) be linear in both time and memory in order to scale to the size of present-day knowledge graphs, and keep up with their growth.

A natural way to handle any possible set of relations is to use the classic canonical polyadic (CP) decomposition Hitchcock (1927), which yields two different embeddings for each entity and thus low prediction performances as shown in Section 5. With unique entity embeddings, multilinear products scale well and can naturally handle both symmetry and (ir)-reflexivity of relations, and when combined with an appropriate loss function, dot products can even handle transitivity Bouchard et al. (2015). However, dealing with antisymmetric—and more generally asymmetric—relations has so far almost always implied superlinear time and space complexity Nickel et al. (2011); Socher et al. (2013) (see Section 2), making models prone to overfitting and not scalable. Finding the best trade-off between expressiveness, generalization and complexity is the keystone of embedding models.

This paper extends a previously published article Trouillon et al. (2016). This extended version adds proofs of existence of the proposed model in both single and multi-relational settings, as well as proofs of the non-uniqueness of the complex embeddings for a given relation. Bounds on the rank of the proposed decomposition are also demonstrated and discussed. The learning algorithm is provided in more details, and more experiments are provided, especially regarding the training time of the models.

The remainder of the paper is organized as follows. We first provide justification and intuition for using complex embeddings in the square matrix case (Section 2), where there is only a single type of relation between entities, and show the existence of the proposed decomposition for all possible relations. The formulation is then extended to a stacked set of square matrices in a third-order tensor to represent multiple relations (Section 3). The stochastic gradient descent algorithm used to learn the model is detailed in Section 4, where we present an equivalent reformulation of the proposed model that involves only real embeddings. This should help practitioners when implementing our method, without requiring the use of complex numbers in their software implementation. We then describe experiments on large-scale public benchmark knowledge graphs in which we empirically show that this representation leads not only to simpler and faster algorithms, but also gives a systematic accuracy improvement over current state-of-the-art alternatives (Section 5). Related work is discussed in Section 6.

Relations as the Real Parts of Low-Rank Normal Matrices

We consider in this section a simplified link prediction task with a single relation, and introduce complex embeddings for low-rank matrix factorization.

We will first discuss the desired properties of embedding models, show how this problem relates to the spectral theorems, and discuss the classes of matrices these theorems encompass in the real and in the complex case. We then propose a new matrix decomposition—to the best of our knowledge—and a proof of its existence for all real square matrices. Finally we discuss the rank of the proposed decomposition.

Let $\mathcal{E}$ be a set of entities, with $|\mathcal{E}|=n$. The truth of the single relation holding between two entities is represented by a sign value $y_{so}\in\{-1,1\}$, where 1 represents true facts and -1 false facts, $s\in\mathcal{E}$ is the subject entity and $o\in\mathcal{E}$ is the object entity. The probability for the relation holding true is given by

In this work we pursue three objectives: finding a generic structure for $X$ that leads to $(i)$ a computationally efficient model, $(ii)$ an expressive enough approximation of common relations in real world knowledge graphs, and $(iii)$ good generalization performances in practice. Standard matrix factorization approximates $X$ by a matrix product $UV^{\top}$, where $U$ and $V$ are two functionally-independent $n\times K$ matrices, $K$ being the rank of the matrix. Within this formulation it is assumed that entities appearing as subjects are different from entities appearing as objects. In the Netflix challenge Koren et al. (2009) for example, each row $u_{i}$ corresponds to the user $i$ and each column $v_{j}$ corresponds to the movie $j$. This extensively studied type of model is closely related to the singular value decomposition (SVD) and fits well with the case where the matrix $X$ is rectangular.

However, in many knowledge graph completion problems, the same entity $i$ can appear as both subject or object and will have two different embedding vectors, $u_{i}$ and $v_{i}$, depending on whether it appears as subject or object of a relation. It seems natural to learn unique embeddings of entities, as initially proposed by Nickel et al. (2011) and Bordes et al. (2011) and since then used systematically in other prominent approaches Bordes et al. (2013b); Yang et al. (2015); Socher et al. (2013). In the factorization setting, using the same embeddings for left- and right-side factors boils down to a specific case of eigenvalue decomposition: orthogonal diagonalization.

The spectral theorem for symmetric matrices tells us that a matrix is orthogonally diagonalizable if and only if it is symmetric Cauchy (1829). It is therefore often used to approximate covariance matrices, kernel functions and distance or similarity matrices.

However as previously stated, this paper is explicitly interested in problems where matrices—and thus the relation patterns they represent—can also be antisymmetric, or even not have any particular symmetry pattern at all (asymmetry). In order to both use a unique embedding for entities and extend the expressiveness to asymmetric relations, researchers have generalised the notion of dot products to scoring functions, also known as composition functions, that allow more general combinations of embeddings. We briefly recall several examples of scoring functions in Table 2.1.1, as well as the extension proposed in this paper.

These models propose different trade-offs between the three essential points:

Expressiveness, which is the ability to represent symmetric, antisymmetric and more generally asymmetric relations.

Scalability, which means keeping linear time and space complexity scoring function.

Generalization, for which having unique entity embeddings is critical.

RESCAL Nickel et al. (2011) and NTN Socher et al. (2013) are very expressive, but their scoring functions have quadratic complexity in the rank of the factorization. More recently the HolE model Nickel et al. (2016b) proposes a solution that has quasi-linear complexity in time and linear space complexity. DistMult Yang et al. (2015) can be seen as a joint orthogonal diagonalization with real embeddings, hence handling only symmetric relations. Conversely, TransE Bordes et al. (2013b) handles symmetric relations to the price of strong constraints on its embeddings. The canonical-polyadic decomposition (CP) Hitchcock (1927) generalizes poorly with its different embeddings for entities as subject and as object.

We reconcile expressiveness, scalability and generalization by going back to the realm of well-studied matrix factorizations, and making use of complex linear algebra, a scarcely used tool in the machine learning community.

In general, there are many other possible couples of matrices $E$ and $W$ that preserve the real part of the decomposition. In practice however this is no synonym of low generalization abilities, as many effective matrix and tensor decomposition methods used in machine learning lead to non-unique solutions Paatero and Tapper (1994); Nickel et al. (2011). In this case also, the learned representations prove useful as shown in the experimental section.

Addressing knowledge graph completion with data-driven approaches assumes that there is a sufficient regularity in the observed data to generalize to unobserved facts. When formulated as a matrix completion problem, as it is the case in this section, one way of implementing this hypothesis is to make the assumption that the matrix has a low rank or approximately low rank. We first discuss the rank of the proposed decomposition, and then introduce the sign-rank and extend the bound developed on the rank to the sign-rank.

First, we recall one definition of the rank of a matrix Horn and Johnson (2012).

We now show that any rank $k$ real square matrix can be reconstructed from a $2k$-dimensional unitary diagonalization.

Since $E$ is not necessarily square, we replace the unitary requirement of Corollary 1 by the requirement that its columns form an orthonormal basis of its smallest dimension, $2k$.

2.2 Sign-Rank Upper Bound

Since we encode the truth values of each fact with $\pm 1$, we deal with square sign matrices: $Y\in\{-1,1\}^{n\times n}$. Sign matrices have an alternative rank definition, the sign-rank.

where the value $c=0$ is here arbitrarily assigned to $1$ to allow zero entries in $X$, conversely to the stricter usual definition of the sign-rank.

To make generalization possible, we hypothesize that the true matrix $Y$ has a low sign-rank, and thus can be reconstructed by the sign of a low-rank score matrix $X$. The low sign-rank assumption is theoretically justified by the fact that the sign-rank is a natural complexity measure of sign matrices Linial et al. (2007) and is linked to learnability Alon et al. (2016) and empirically confirmed by the wide success of factorization models Nickel et al. (2016a).

Using Corollary 2, we can now show that any square sign matrix of sign-rank $k$ can be reconstructed from a rank $2k$ unitary diagonalization.

Previous attempts to approximate the sign-rank in relational learning did not use complex numbers. They showed the existence of compact factorizations under conditions on the sign matrix Nickel et al. (2014), or only in specific cases Bouchard et al. (2015). In contrast, our results show that if a square sign matrix has sign-rank $k$, then it can be exactly decomposed through a $2k$-dimensional unitary diagonalization.

Finding the $K$ that matches the sign-rank of $Y$ corresponds to finding the smallest $K$ that brings the 0–1 loss on $X$ to , as link prediction can be seen as binary classification of the facts. In practice, and as classically done in machine learning to avoid this NP-hard problem, we use a continuous surrogate of the 0–1 loss, in this case the logistic loss as described in Section 4, and validate models on different values of $K$, as described in Section 5.

2.3 Rank Bound Discussion

Corollaries 2 and 3 use the aforementioned subadditive property of the rank to derive the $2k$ upper bound. Let us give an example for which this bound is strictly greater than $k$.

Consider the following $2$-by-$2$ sign matrix:

Not only is this matrix not normal, but one can also easily check that there is no real normal $2$-by-$2$ matrix that has the same sign-pattern as $Y$. Clearly, $Y$ is a rank $1$ matrix since its columns are linearly dependent, hence its sign-rank is also $1$. From Corollary 3, we know that there is a normal matrix whose real part has the same sign-pattern as $Y$, and whose rank is at most $2$.

This example shows that the upper bound on the rank of the unitary diagonalization showed in Corollaries 2 and 3 can be strictly greater than $k$, the rank or sign-rank, of the decomposed matrix. However, there might be other examples for which the addition of an imaginary part could—additionally to making the matrix normal—create some linear dependence between the rows/columns and thus decrease the rank of the matrix, up to a factor of 2.

We summarize this section in three points:

The proposed factorization encompasses all possible score matrices $X$ for a single binary relation.

By construction, the factorization is well suited to represent both symmetric and antisymmetric relations.

Relation patterns can be efficiently approximated with a low-rank factorization using complex-valued embeddings.

Given one relation $r\in\mathcal{R}$ and two entities $s,o$ $\in\mathcal{E}$, the probability that the fact r(s,o) is true given by:

where $\phi$ is the scoring function of the model considered and $\Theta$ denotes the model parameters. We denote the set of all possible facts (or triples) for a knowledge graph by $\mathcal{T}=\mathcal{R}\times\mathcal{E}\times\mathcal{E}$. While the tensor $\mathbf{X}$ as a whole is unknown, we assume that we observe a set of true and false triples $\Omega=\{((r,s,o),y_{rso})\,|\,(r,s,o)\in\mathcal{T}_{\Omega}\}$ where $y_{rso}\in\{-1,1\}$ and $\mathcal{T}_{\Omega}\subseteq\mathcal{T}$ is the set of observed triples. The goal is to find the probabilities of entries $y_{r^{\prime}s^{\prime}o^{\prime}}$ for a set of targeted unobserved triples $\{(r^{\prime},s^{\prime},o^{\prime})\in\mathcal{T}\setminus\mathcal{T}_{\Omega}\}$.

Depending on the scoring function $\phi(r,s,o;\Theta)$ used to model the score tensor $\mathbf{X}$, we obtain different models. Examples of scoring functions are given in Table 2.1.1.

These equations provide two interesting views of the model:

Changing the representation: Equation (6) would correspond to DistMult with real embeddings (see Table 2.1.1), but handles asymmetry thanks to the complex conjugate of the object-entity embedding.

Changing the scoring function: Equation (7) only involves real vectors corresponding to the real and imaginary parts of the embeddings and relations.

From a geometrical point of view, each relation embedding $w_{r}$ is an anisotropic scaling of the basis defined by the entity embeddings $E$, followed by a projection onto the real subspace.

2 Existence of the Tensor Factorization

Let us first discuss the existence of the multi-relational model where the rank of the decomposition $K\leq n$, which relates to simultaneous unitary decomposition.

To apply Theorem 3 to the proposed factorization, we would have to make the hypothesis that the relation score matrices $X_{r}$ are a commuting family, which is too strong a hypothesis. Actually, the model is slightly different since we take only the real part of the tensor factorization. In the single-relation case, taking only the real part of the decomposition rids us of the normality requirement of Theorem 1 for the decomposition to exist, as shown in Theorem 2.

In the multiple-relation case, it is an open question whether taking the real part of the simultaneous unitary diagonalization will enable us to decompose families of arbitrary real square matrices—that is with a single unitary matrix $E$ that has at most $n$ columns. Though it seems unlikely, we could not find a counter-example yet.

However, by letting the rank of the tensor factorization $K$ to be greater than $n$, we can show that the proposed tensor decomposition exists for families of arbitrary real square matrices, by simply concatenating the decomposition of Theorem 2 of each real square matrix $X_{i}$.

Let $E=\left[E_{1}\ldots E_{m}\right]$, and

where $\mathbf{0}^{l\times l}$ the zero $l\times l$ matrix. Therefore $X_{i}=(E\Lambda_{i}E^{*})$ for all $i$ in $1,\ldots,m$.

By construction, the rank of the decomposition is at most $nm$. When $m\leq n$, this bound actually matches the general upper bound on the rank of the canonical polyadic (CP) decomposition Hitchcock (1927); Kruskal (1989). Since $m$ corresponds to the number of relations and $n$ to the number of entities, $m$ is always smaller than $n$ in real world knowledge graphs, hence the bound holds in practice.

Though when it comes to relational learning, we might expect the actual rank to be much lower than $nm$ for two reasons. The first one, as discussed above, is that we are dealing with sign tensors, hence the rank of the matrices $X_{r}$ need only match the sign-rank of the partially-observed matrices $Y_{r}$. The second one is that the matrices are related to each other, as they all represent the same entities in different relations, and thus benefit from sharing latent dimensions. As opposed to the construction exposed in the proof of Theorem 4, where other relations dimensions are canceled out. In practice, the rank needed to generalize well is indeed much lower than $nm$ as we show experimentally in Figure 5.

Also, note that with the construction of the proof of Theorem 4, the matrix $E=\left[E_{1}\ldots E_{m}\right]$ is not unitary any more. However the unitary constraints in the matrix case serve only the proof of existence, which is just one solution among the infinite ones of same rank. In practice, imposing orthonormality is essentially a numerical commodity for the decomposition of dense matrices, through iterative methods for example Saad (1992). When it comes to matrix and tensor completion, and thus generalisation, imposing such constraints is more of a numerical hassle than anything else, especially for gradient methods. As there is no apparent link between orthonormality and generalisation properties, we did not impose these constraints when learning this model in the following experiments.

Algorithm 1 describes stochastic gradient descent (SGD) to learn the proposed multi-relational model with the AdaGrad learning-rate updates Duchi et al. (2011). We refer to the proposed model as ComplEx, for Complex Embeddings. We expose a version of the algorithm that uses only real-valued vectors, in order to facilitate its implementation. To do so, we use separate real-valued representations of the real and imaginary parts of the embeddings.

These real and imaginary part vectors are initialized with vectors having a zero-mean normal distribution with unit variance. If the training set $\Omega$ contains only positive triples, negatives are generated for each batch using the local closed-world assumption as in Bordes et al. (2013b). That is, for each triple, we randomly change either the subject or the object, to form a negative example. In this case the parameter $\eta>0$ sets the number of negative triples to generate for each positive triple. Collision with positive triples in $\Omega$ is not checked, as it occurs rarely in real world knowledge graphs as they are largely sparse, and may also be computationally expensive.

Squared gradients are accumulated to compute AdaGrad learning rates, then gradients are updated. Every $s$ iterations, the parameters $\Theta$ are evaluated over the evaluation set $\Omega_{v}$ (evaluate_AP_or_MRR$(\Omega_{v};\Theta)$ function in Algorithm 1). If the data set contains both positive and negative examples, average precision (AP) is used to evaluate the model. If the data set contains only positives, then mean reciprocal rank (MRR) is used as average precision cannot be computed without true negatives. The optimization process is stopped when the measure considered decreases compared to the last evaluation (early stopping).

Bern($p$) is the Bernoulli distribution, the one_random_sample$(\mathcal{E})$ function sample uniformly one entity in the set of all entities $\mathcal{E}$, and the sample_batch_of_size_b$(\Omega,b)$ function sample $b$ true and false triples uniformly at random from the training set $\Omega$.

where each entity and each relation has two real embeddings.

where $\odot$ is the element-wise (Hadamard) product.

We optimized the negative log-likelihood of the logistic model described in Equation (5) with $L^{2}$ regularization on the parameters $\Theta$:

To handle regularization, note that using separate representations for the real and imaginary parts does not change anything as the squared $L^{2}$-norm of a complex vector $v=v^{\prime}+iv^{\prime\prime}$ is the sum of the squared modulus of each entry:

which is actually the sum of the $L^{2}$-norms of the vectors of the real and imaginary parts.

We can finally write the gradient of $\gamma$ with respect to a real embedding $v$ for one triple $(r,s,o)$ and its truth value $y$:

We evaluated the method proposed in this paper on both synthetic and real data sets. The synthetic data set contains both symmetric and antisymmetric relations, whereas the real data sets are standard link prediction benchmarks based on real knowledge graphs.

We compared ComplEx to state-of-the-art models, namely TransE Bordes et al. (2013b), DistMult Yang et al. (2015), RESCAL Nickel et al. (2011) and also to the canonical polyadic decomposition (CP) Hitchcock (1927), to emphasize empirically the importance of learning unique embeddings for entities. For experimental fairness, we reimplemented these models within the same framework as the ComplEx model, using a Theano-based SGD implementationhttps://github.com/lmjohns3/downhill Bergstra et al. (2010).

For the TransE model, results were obtained with its original max-margin loss, as it turned out to yield better results for this model only. To use this max-margin loss on data sets with observed negatives (Sections 5.1 and 5.2), positive triples were replicated when necessary to match the number of negative triples, as described in Garcia-Duran et al. (2016). All other models are trained with the negative log-likelihood of the logistic model (Equation (10)). In all the following experiments we used a maximum number of iterations $m=1000$, a batch size $b=\frac{|\Omega|}{100}$, and validated the models for early stopping every $s=50$ iterations.

To assess our claim that ComplEx can accurately model jointly symmetry and antisymmetry, we randomly generated a knowledge graph of two relations and 30 entities. One relation is entirely symmetric, while the other is completely antisymmetric. This data set corresponds to a $2\times 30\times 30$ tensor. Figure 2 shows a part of this randomly generated tensor, with a symmetric slice and an antisymmetric slice, decomposed into training, validation and test sets. To ensure that all test values are predictable, the upper triangular parts of the matrices are always kept in the training set, and the diagonals are unobserved. We conducted a 5-fold cross-validation on the lower-triangular matrices, using the upper-triangular parts plus 3 folds for training, one fold for validation and one fold for testing. Each training set contains 1392 observed triples, whereas validation and test sets contain 174 triples each.

Figure 3 shows the best cross-validated average precision (area under the precision-recall curve) for different factorization models of ranks ranging up to 50. The regularization parameter $\lambda$ is validated in $\{$0.1, 0.03, 0.01, 0.003,0.001, 0.0003, 0.00001, 0.0$\}$ and the learning rate $\alpha$ was initialized to 0.5.

As expected, DistMult Yang et al. (2015) is not able to model antisymmetry and only predicts the symmetric relations correctly. Although TransE Bordes et al. (2013b) is not a symmetric model, it performs poorly in practice, particularly on the antisymmetric relation. RESCAL Nickel et al. (2011), with its large number of parameters, quickly overfits as the rank grows. Canonical Polyadic (CP) decomposition Hitchcock (1927) fails on both relations as it has to push symmetric and antisymmetric patterns through the entity embeddings. Surprisingly, only ComplEx succeeds even on such simple data.

2 Real Fully-Observed Data Sets: Kinships and UMLS

We then compare all models on two fully observed data sets, that contain both positive and negative triples, also called the closed-world assumption. The Kinships data set Denham (1973) describes the 26 different kinship relations of the Alyawarra tribe in Australia, among 104 individuals. The unified medical language system (UMLS) data set McCray (2003) represents 135 medical concepts and diseases, linked by 49 relations describing their interactions. Metadata for the two data sets is summarized in Table 2.

We performed a 10-fold cross-validation, keeping 8 for training, one for validation and one for testing. Figure 4 shows the best cross-validated average precision for ranks ranging up to 50, and error bars show the standard deviation over the 10 runs. The regularization parameter $\lambda$ is validated in $\{$0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.00001, 0.0$\}$ and the learning rate $\alpha$ was initialized to 0.5.

On both data sets ComplEx, RESCAL and CP are very close, with a slight advantage for ComplEx on Kinships, and for RESCAL on UMLS. DistMult performs poorly here as many relations are antisymmetric both in UMLS (causal links, anatomical hierarchies) and Kinships (being father, uncle or grand-father).

The fact that CP, RESCAL and ComplEx work so well on these data sets illustrates the importance of having an expressive enough model, as DistMult fails because of antisymmetry; the power of the multilinear product—that is the tensor factorization approach—as TransE can be seen as a sum of bilinear products Garcia-Duran et al. (2016); but not yet the importance of having unique entity embeddings, as CP works well. We believe having separate subject and object-entity embeddings works well under the closed-world assumption, because of the amount of training data compared to the number of embeddings to learn. Though when only a fractions of the positive training examples are observed (as it is most often the case), we will see in the next experiments that enforcing unique entity embeddings is key to good generalization.

3 Real Sparse Data Sets: FB15K and WN18

Finally, we evaluated ComplEx on the FB15K and WN18 data sets, as they are well established benchmarks for the link prediction task. FB15K is a subset of Freebase Bollacker et al. (2008), a curated knowledge graph of general facts, whereas WN18 is a subset of WordNet Fellbaum (1998), a database featuring lexical relations between words. We used the same training, validation and test set splits as in Bordes et al. (2013b). Table 3 summarizes the metadata of the two data sets.

As both data sets contain only positive triples, we generated negative samples using the local closed-world assumption, as described in Section 4. For evaluation, we measure the quality of the ranking of each test triple among all possible subject and object substitutions : $r(s^{\prime},o)$ and $r(s,o^{\prime})$, for each $s^{\prime},o^{\prime}$ in $\mathcal{E}$, as used in previous studies Bordes et al. (2013b); Nickel et al. (2016b). Mean Reciprocal Rank (MRR) and Hits at $N$ are standard evaluation measures for these data sets and come in two flavours: raw and filtered. The filtered metrics are computed after removing all the other positive observed triples that appear in either training, validation or test set from the ranking, whereas the raw metrics do not remove these.

Since ranking measures are used, previous studies generally preferred a max-margin ranking loss for the task Bordes et al. (2013b); Nickel et al. (2016b). We chose to use the negative log-likelihood of the logistic model—as described in the previous section—as it is a continuous surrogate of the sign-rank, and has been shown to learn compact representations for several important relations, especially for transitive relations Bouchard et al. (2015). As previously stated, we tried both losses in preliminary work, and indeed training the models with the log-likelihood yielded better results than with the max-margin ranking loss, especially on FB15K—except with TransE.

We report both filtered and raw MRR, and filtered Hits at 1, 3 and 10 in Table 4 for the evaluated models. The HolE model has recently been shown to be equivalent to ComplEx Hayashi and Shimbo (2017), we record the original results for HolE as reported in Nickel et al. (2016b) and briefly discuss the discrepancy of results obtained with ComplEx.

Reported results are given for the best set of hyper-parameters evaluated on the validation set for each model, after a distributed grid-search on the following values: $K\in\{$10, 20, 50, 100, 150, 200$\}$, $\lambda\in\{$0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.0$\}$, $\alpha\in\{$1.0, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01$\}$, $\eta\in\{$1, 2, 5, 10$\}$ with $\lambda$ the $L^{2}$ regularization parameter, $\alpha$ the initial learning rate, and $\eta$ the number of negatives generated per positive training triple. We also tried varying the batch size but this had no impact and we settled with 100 batches per epoch. With the best hyper-parameters, training the ComplEx model on a single GPU (NVIDIA Tesla P40) takes 45 minutes on WN18 ($K=150,\eta=1$), and three hours on FB15K ($K=200,\eta=10$).

3.2 Results

WN18 describes lexical and semantic hierarchies between concepts and contains many antisymmetric relations such as hypernymy, hyponymy, and being part of. Indeed, the DistMult and TransE models are outperformed here by ComplEx and HolE, which are on a par with respective filtered MRR scores of 0.941 and 0.938, which is expected as both models are equivalent.

Table 5 shows the filtered MRR for the reimplemented models and each relation of WN18, confirming the advantage of ComplEx on antisymmetric relations while losing nothing on the others. 2D projections of the relation embeddings (Figures 8 & 9) visually corroborate the results.

On FB15K, the gap is much more pronounced and the ComplEx model largely outperforms HolE, with a filtered MRR of 0.692 and 59.9% of Hits at 1, compared to 0.524 and 40.2% for HolE. This difference of scores between the two models, though they have been proved to be equivalent Hayashi and Shimbo (2017), is due to the use of the aforementioned max-margin loss in the original HolE publication Nickel et al. (2016b) that performs worse than the log-likelihood on this dataset, and to the generation of more than one negative sample per positive in these experiments. This has been confirmed and discussed in details by Trouillon and Nickel (2017). The fact that DistMult yields fairly high scores (0.654 filtered MRR) is also due to the task itself and the evaluation measures used. As the dataset only involves true facts, the test set never includes the opposite facts $r(o,s)$ of each test fact $r(s,o)$ for antisymmetric relations—as the opposite fact is always false. Thus highly scoring the opposite fact barely impacts the rankings for antisymmetric relations. This is not the case in the fully observed experiments (Section 5.2), as the opposite fact is known to be false—for antisymmetric relations—and largely impacts the average precision of the DistMult model (Figure 4).

RESCAL, that represents each relation with a $K\times K$ matrix, performs well on WN18 as there are few relations and hence not so many parameters. On FB15K though, it probably overfits due to the large number of relations and thus the large number of parameters to learn, and performs worse than a less expressive model like DistMult. On both data sets, TransE and CP are largely left behind. This illustrates again the power of the multilinear product in the first case, and the importance of learning unique entity embeddings in the second. CP performs especially poorly on WN18 due to the small number of relations, which magnifies this subject/object difference.

Figure 5 shows that the filtered MRR of the ComplEx model quickly converges on both data sets, showing that the low-rank hypothesis is reasonable in practice. The little gain of performances for ranks comprised between $50$ and $200$ also shows that ComplEx does not perform better because it has twice as many parameters for the same rank—the real and imaginary parts—compared to other linear space complexity models but indeed thanks to its better expressiveness.

Best ranks were generally 150 or 200, in both cases scores were always very close for all models, suggesting there was no need to grid-search on higher ranks. The number of negative samples per positive sample also had a large influence on the filtered MRR on FB15K (up to +0.08 improvement from 1 to 10 negatives), but not much on WN18. On both data sets regularization was important (up to +0.05 on filtered MRR between $\lambda=0$ and optimal one). We found the initial learning rate to be very important on FB15K, while not so much on WN18. We think this may also explain the large gap of improvement ComplEx provides on this data set compared to previously published results—as DistMult results are also better than those previously reported Yang et al. (2015)—along with the use of the log-likelihood objective. It seems that in general AdaGrad is relatively insensitive to the initial learning rate, perhaps causing some overconfidence in its ability to tune the step size online and consequently leading to less efforts when selecting the initial step size.

4 Training time

As defended in Section 2, having a linear time and space complexity becomes critical when the dataset grows. To illustrate this, we report in Figure 6 the evolution of the filtered MRR on the validation set as a function of time, for the best set of validated hyper-parameters for each model. The convergence criterion used is the decrease of the validation filtered MRR—computed every 50 iterations—with a maximum number of iterations of 1000 (see Algorithm 1). All models have a linear complexity except for RESCAL that has a quadratic one in the rank of the decomposition, as it learns one matrix embedding for each relation $r\in\mathcal{R}$. Timings are measured on a single NVIDIA Tesla P40 GPU.

On WN18, all models reach convergence in a reasonable time, between 15 minutes and 1 hour and 20 minutes. The difference between RESCAL and the other models is not sharp there, first because its optimal embedding size ($K=50$) is lower compared to the other models. Secondly, there are only $|\mathcal{R}|=18$ relations in WN18, hence the memory footprint of RESCAL is pretty similar to the other models—because it represents only relations with matrices and not entities. On FB15K, the difference is much more pronounced, as RESCAL optimal rank is similar to the other models; and with $|\mathcal{R}|=1345$ relations, RESCAL has a much higher memory footprint, which implies more processor cache misses due to the uniformly-random nature of the SGD sampling.

RESCAL took more than four days to train on FB15K, whereas other models took between 40 minutes and 3 hours. While a few days might seem manageable, this could not be the case on larger data sets, as FB15K is but a small subset of Freebase that contains $|\mathcal{R}|=35000$ relations Bollacker et al. (2008). This experimentally supports our claim that linear complexity is required for scalability.

We further investigated the influence of the number of negatives generated per positive training sample. In the previous experiment, due to computational limitations, the number of negatives per training sample, $\eta$, was validated over the set $\{1,2,5,10\}$. On WN18 it proved to be of no help to have more than one generated negative per positive. Here we explore in which proportions increasing the number of generated negatives leads to better results on FB15K. To do so, we fixed the best validated $\lambda,K,\alpha$ obtained from the previous experiment. We then let $\eta$ vary in $\{1,2,5,10,20,50,100,200\}$.

Figure 7 shows the influence of the number of generated negatives per positive training triple on the performance of ComplEx on FB15K. Generating more negatives clearly improves the results up to 100 negative triples, with a filtered MRR of 0.737 and 64.8% of Hits@1, before decreasing again with 200 negatives, probably due to the too large class imbalance. The model also converges with fewer epochs, which compensates partially for the additional training time per epoch, up to 50 negatives. It then grows linearly as the number of negatives increases.

4.2 WN18 Embeddings Visualization

We used principal component analysis (PCA) to visualize embeddings of the relations of the WordNet data set (WN18). We plotted the four first components of the best DistMult and ComplEx model’s embeddings in Figures 8 & 9. For the ComplEx model, we simply concatenated the real and imaginary parts of each embedding.

Most of WN18 relations describe hierarchies, and are thus antisymmetric. Each of these hierarchic relations has its inverse relation in the data set. For example: hypernym / hyponym, part_of / has_part, synset_domain_topic_of / member_of_domain_topic. Since DistMult is unable to model antisymmetry, it will correctly represent the nature of each pair of opposite relations, but not the direction of the relations. Loosely speaking, in the hypernym / hyponym pair the nature is sharing semantics, and the direction is that one entity generalizes the semantics of the other. This makes DistMult representing the opposite relations with very close embeddings. It is especially striking for the third and fourth principal component (Figure 9). Conversely, ComplEx manages to oppose spatially the opposite relations.

We first discuss related work about complex-valued matrix and tensor decompositions, and then review other approaches for knowledge graph completion.

When factorization methods are applied, the representation of the decomposition is generally chosen in accordance with the data, despite the fact that most real square matrices only have eigenvalues in the complex domain. Indeed in the machine learning community, the data is usually real-valued, and thus eigendecomposition is used for symmetric matrices, or other decompositions such as (real-valued) singular value decomposition Beltrami (1873), non-negative matrix factorization Paatero and Tapper (1994), or canonical polyadic decomposition when it comes to tensors Hitchcock (1927).

Conversely, in signal processing, data is often complex-valued Stoica and Moses (2005) and the complex-valued counterparts of these decompositions are then used. Joint diagonalization is also a much more common tool than in machine learning for decomposing sets of (complex) dense square matrices Belouchrani et al. (1997); De Lathauwer et al. (2001).

Some works on recommender systems use complex numbers as an encoding facility, to merge two real-valued relations, similarity and liking, into one single complex-valued matrix which is then decomposed with complex embeddings Kunegis et al. (2012); Xie et al. (2015). Still, unlike our work, it is not real data that is decomposed in the complex domain.

In deep learning, Danihelka et al. (2016) proposed an LSTM extended with an associative memory based on complex-valued vectors for memorization tasks, and Hu et al. (2016) a complex-valued neural network for speech synthesis. In both cases again, the data is first encoded in complex vectors that are then fed into the network.

Conversely to these contributions, this work suggests that processing real-valued data with complex-valued representation, through a projection onto the real-valued subspace, can be a very simple way of increasing the expressiveness of the model considered.

2 Knowledge Graph Completion

Many knowledge graphs have recently arisen, pushed by the W3C recommendation to use the resource description framework (RDF) Cyganiak et al. (2014) for data representation. Examples of such knowledge graphs include DBPedia Auer et al. (2007), Freebase Bollacker et al. (2008) and the Google Knowledge Vault Dong et al. (2014). Motivating applications of knowledge graph completion include question answering Bordes et al. (2014b) and more generally probabilistic querying of knowledge bases Huang and Liu (2009); Krompaß et al. (2014).

First approaches to relational learning relied upon probabilistic graphical models Getoor and Taskar (2007), such as bayesian networks Friedman et al. (1999) and markov logic networks Richardson and Domingos (2006); Raedt et al. (2016).

With the first embedding models, asymmetry of relations was quickly seen as a problem and asymmetric extensions of tensors were studied, mostly by either considering independent embeddings Franz et al. (2009) or considering relations as matrices instead of vectors in the RESCAL model Nickel et al. (2011), or both Sutskever (2009). Direct extensions were based on uni-,bi- and trigram latent factors for triple data Garcia-Duran et al. (2016), as well as a low-rank relation matrix Jenatton et al. (2012). Bordes et al. (2014a) propose a two-layer model where subject and object embeddings are first separately combined with the relation embedding, then each intermediate representation is combined into the final score.

Pairwise interaction models were also considered to improve prediction performances. For example, the Universal Schema approach Riedel et al. (2013) factorizes a 2D unfolding of the tensor (a matrix of entity pairs vs. relations) while Welbl et al. (2016) extend this also to other pairs. Riedel et al. (2013) also consider augmenting the knowledge graph facts by exctracting them from textual data, as does Toutanova et al. (2015). Injecting prior knowledge in the form of Horn clauses in the objective loss of the Universal Schema model has also been considered Rocktaschel et al. (2015). Chang et al. (2014) enhance the RESCAL model to take into account information about the entity types. For recommender systems (thus with different subject/object sets of entities), Baruch (2014) proposed a non-commutative extension of the CP decomposition model. More recently, Gaifman models that learn neighborhood embeddings of local structures in the knowledge graph showed competitive performances Niepert (2016).

In the Neural Tensor Network (NTN) model, Socher et al. (2013) combine linear transformations and multiple bilinear forms of subject and object embeddings to jointly feed them into a nonlinear neural layer. Its non-linearity and multiple ways of including interactions between embeddings gives it an advantage in expressiveness over models with simpler scoring function like DistMult or RESCAL. As a downside, its very large number of parameters can make the NTN model harder to train and overfit more easily.

The original multilinear DistMult model is symmetric in subject and object for every relation Yang et al. (2015) and achieves good performance on FB15K and WN18 data sets. However it is likely due to the absence of true negatives in these data sets, as discussed in Section 5.3.2.

The TransE model from Bordes et al. (2013b) also embeds entities and relations in the same space and imposes a geometrical structural bias into the model: the subject entity vector should be close to the object entity vector once translated by the relation vector.

A recent novel way to handle antisymmetry is via the Holographic Embeddings (HolE) model by Nickel et al. (2016b). In HolE the circular correlation is used for combining entity embeddings, measuring the covariance between embeddings at different dimension shifts. This model has been shown to be equivalent to the ComplEx model Hayashi and Shimbo (2017); Trouillon and Nickel (2017).

Though the decomposition proposed in this paper is clearly not unique, it is able to learn meaningful representations. Still, characterizing all possible unitary diagonalizations that preserve the real part is an interesting open question. Especially in an approximation setting with a constrained rank, in order to characterize the decompositions that minimize a given reconstruction error. That might allow the creation of an iterative algorithm similar to eigendecomposition iterative methods Saad (1992) for computing such a decomposition for any given real square matrix.

The proposed decomposition could also find applications in many other asymmetric square matrices decompositions applications, such as spectral graph theory for directed graphs Cvetković et al. (1997), but also factorization of asymmetric measures matrices such as asymmetric distance matrices Mao and Saul (2004) and asymmetric similarity matrices Pirasteh et al. (2015).

From an optimization point of view, the objective function (Equation (10)) is clearly non-convex, and we could indeed not be reaching a globally optimal decomposition using stochastic gradient descent. Recent results show that there are no spurious local minima in the completion problem of positive semi-definite matrix Ge et al. (2016); Bhojanapalli et al. (2016). Studying the extensibility of these results to our decomposition is another possible line of future work. The first step would be generalizing these results to symmetric real-valued matrix completion, then generalization to normal matrices should be straightforward. The two last steps would be extending to matrices that are expressed as real part of normal matrices, and finally to the joint decomposition of such matrices as a tensor. We indeed noticed a remarkable stability of the scores across different random initialization of ComplEx for the same hyper-parameters, which suggests the possibility of such theoretical property.

Practically, an obvious extension is to merge our approach with known extensions to tensor factorization models in order to further improve predictive performance. For example, the use of pairwise embeddings Riedel et al. (2013); Welbl et al. (2016) together with complex numbers might lead to improved results in many situations that involve non-compositionality. Adding bigram embeddings to the objective could also improve the results as shown on other models Garcia-Duran et al. (2016). Another direction would be to develop a more intelligent negative sampling procedure, to generate more informative negatives with respect to the positive sample from which they have been sampled. This would reduce the number of negatives required to reach good performance, thus accelerating training time. Extension to relations between more than two entities, $n$-tuples, is not straightforward, as ComplEx’s expressiveness comes from the complex conjugation of the object-entity, that breaks the symmetry between the subject and object embeddings in the scoring function. This stems from the Hermitian product, which seems to have no standard multilinear extension in the linear algebra literature, this question hence remains largely open.

We described a new matrix and tensor decomposition with complex-valued latent factors called ComplEx. The decomposition exists for all real square matrices, expressed as the real part of normal matrices. The result extends to sets of real square matrices—tensors—and answers to the requirements of the knowledge graph completion task : handling a large variety of different relations including antisymmetric and asymmetric ones, while being scalable. Experiments confirm its theoretical versatility, as it substantially improves over the state-of-the-art on real knowledge graphs. It shows that real world relations can be efficiently approximated as the real part of low-rank normal matrices. The generality of the theoretical results and the effectiveness of the experimental ones motivate for the application to other real square matrices factorization problems. More generally, we hope that this paper will stimulate the use of complex linear algebra in the machine learning community, even and especially for processing real-valued data.

This work was supported in part by the Association Nationale de la Recherche et de la Technologie through the CIFRE grant 2014/0121, in part by the Paul Allen Foundation through an Allen Distinguished Investigator grant, and in part by a Google Focused Research Award. We would like to thank Ariadna Quattoni, Stéphane Clinchant, Jean-Marc Andréoli, Sofia Michel, Alejandro Blumentals, Léo Hubert and Pierre Comon for their helpful comments and feedback.