{'fraction_non_alphanumeric': 0.03729071537290715, 'fraction_numerical': 0.013972602739726028, 'mean_word_length': 4.451543312313309, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Probit Tensor Factorization Model For Relational Learning

With the proliferation of knowledge graphs, modeling data with complex multirelational structure has gained increasing attention in the area of statistical relational learning. One of the most important goals of statistical relational learning is link prediction, i.e., predicting whether certain relations exist in the knowledge graph. A large number of models and algorithms have been proposed to perform link prediction, among which tensor factorization method has proven to achieve state-of-the-art performance in terms of computation efficiency and prediction accuracy. However, a common drawback of the existing tensor factorization models is that the missing relations and non-existing relations are treated in the same way, which results in a loss of information. To address this issue, we propose a binary tensor factorization model with probit link, which not only inherits the computation efficiency from the classic tensor factorization model but also accounts for the binary nature of relational data. Our proposed probit tensor factorization (PTF) model shows advantages in both the prediction accuracy and interpretability


Introduction
With the ever-growing computing technology, a huge volume of structured relational data have been generated in recent years. An example is the creation of knowledge graphs, which are graph-structured knowledge bases that store relations between entities. Notable projects are YAGO (Suchanek et al., 2007), DBpedia (Auer et al., 2007), NELL (Carlson et al., 2010), the Google Knowledge Graph (Singhal, 2012) and so on. These graphs have become an important asset in many aspects. For instance, the Google Knowledge Graph (Singhal, 2012), which currently contains 70 billion relational facts between different entities, provides structured information about how topics are linked to one another.
Such information enables the search engine to provide more relevant results. Knowledge graphs have also been used in several specialized domains. In spite of their huge advantages in representing complex data, these graphs more often than not are very incomplete, noisy and sparse. Consequently, the construction and completion of multi-relational data is of great interest nowadays, where the latter task can be formalized as link prediction in relational learning.
Unlike the majority of research in machine learning that assumes data are independently distributed, this independence assumption does not hold in relational data. In fact, we rely on the dependent relations between data to provide valuable information. The learning model is expected to propagate information about different kinds of relations between a set of entities. Multi-relational data typically describes different types of relations between a set of entities. The entities can be words in natural language processing (NLP) applications, users in social networks, or disease in health record data. One common way to represent the relations between entities is to use triplets of the form (subject, predicate, object). The predicates of interest usually represent the semantic relationship between entities. A formal way to represent the data is to use a three-way binary tensor X of size N × N × K with N as the number of entities and K as the number of relationship types. The k-th slice of the tensor X, which is X k ∈ {−1, 1} N ×N contains information about the k-th relationship, with x ijk = 1, −1 denoting a valid/invalid relationship k from entity i to entity j. The tensor is usually sparse since most of the triplets are unobserved.
Two different assumptions are widely used to interpret the unobserved triplets, including the closed world assumption (CWA) and the open world assumption (OWA). Under the CWA, unobserved triplets are assumed to be invalid, while under the OWA, unobserved triplets can be either valid or invalid. Majority of previous works used the CWA. They did not distinguish between invalid facts and missing facts, which are both encoded as 0. Such a representation will certainly reduce the prediction accuracy because the unobserved facts and invalid facts carry different information. Although most of the large-scale knowledge graphs only contain information for valid triplets, it is not hard to generate invalid samples by using side information. Chang et al. (2014) pointed out that entity type can be used to dictate whether some entities can be legitimate arguments of a given predicate.
Furthermore, most of the previous works on statistical relational learning assumes a Gaussian distribution for the relational data, which gives rise to many variants of the tensor factorization models (Nickel et al., 2011;Chang et al., 2014;Drumond et al., 2016), where the binary nature of the relational data is not taken into account. Such an assumption is neither intuitive nor appropriate from a statistical perspective. Another limitation of these models is that the predictions are continuous scores, which do not have a probabilistic interpretation. A reasonable threshold needs to be manually chosen in order to achieve an accurate prediction.
To modify the Gaussian assumption on the relational data, Nickel and Tresp (2013) proposed a logistic tensor factorization model, RESCAL Logit, where each x ijk 's is assumed to follow a Bernoulli distribution. This model makes the prediction procedure totally data-driven and improves the prediction accuracy. The authors used a quasi-Newton optimization, where some dense matrices need to be calculated in each iteration to get the gradients. This significantly slows down the model fitting and makes it less efficient compared to other tensor factorization models. Besides, Krompaß et al. (2013) extended the RESCAL model with an non-invalid constraint on the parameters.
In this paper, we propose the Probit Tensor Factorization (PTF) Model to address these problems. By imposing a probit link on the binary random variables x ijk 's, we essentially assume that the binary value of each x ijk depends on the sign of a corresponding latent normally distributed random variable Z ijk . The latent Z ijk follows different conditional distributions given x ijk is valid, invalid or missing. All the analyses utilize the R package "ptf", which has been developed for fitting the Probit Tensor Factorization Model. The main contribution of this paper are summarized as follows: First, a PTF Model is proposed to perform link prediction, which not only takes advantage of the computational efficiency from the Gaussian tensor factorization model but also enjoys the prediction accuracy and interpretability of the logistic tensor factorization model.
Second, under the open world assumption (OWA), the unknown triplets are treated differently with known-to-be-invalid triplets, which is a more efficient and proper way of handling missing data. Third, experiments on the moderate-sized real-world datasets demonstrate that our model outperforms the state-of-the-art approaches.
The rest of the paper is organized as follows. First, we summarize several state-of-theart relational learning models. Then we illustrate the framework of our model, with an analysis of computation complexity. After that, we compare prediction performance of our model with the current state-of-art models in real data applications. Finally, we provide a brief conclusion and discussion.

Related Work
In this section, the related work are organized into three parts. We start by reviewing related works in tensor factorization, then summarize related works in the area of relational learning. Finally, we connect our work to the existing network analysis literature.

Related work in tensor factorization
Tensors are generalizations of matrices to higher dimension and have long been studied. It is common to observe missing elements in a multi-dimensional tensor, which motivates the low-rank assumption. The earliest tensors and their decompositions can be traced back to Hitchcock (1927). During the last decades, there has been a large number of researches on tensor factorization related topics targeting on a great variety of applications. For example, Aidini et al. (2018) proposed a 1-bit tensor factorization and applied it to image analysis, Ghadermarzy et al. (2018) generalize the 1-bit tensor factorization model and demonstrated its usage in for context-aware recommender systems.
One of the most widely used tensor factorization is the CANDECOMP/PARAFAC (CP) proposed by Harshman and Lundy (1994), which can be written as where X k is the k-th slice in the tensor whose dimension is I by J; A is an I by R factor loading matrix for Mode A and B is a J by R loading matrix for Mode B; D k is an R by R diagonal matrix; E k is the k-th slice of the residual tensor. Based on the CP decomposition, Rai et al. (2014) proposed a Bayesian low-rank decomposition which is applicable to relational learning setting and achieved good performances.
Another well-known tensor factorization model is the Tucker decomposition (Tucker, 1966). The main idea is to decompose a three-way tensor X ∈ R I×J×K as for i = 1, ..., I, j = 1, ..., J, k = 1, ...K. Xu et al. (2013) proposed a InfTucker models, where binary data and missing data are taken into account.
While these decomposition models have been applied to a great variety of areas (Karami et al., 2012;Rendle and Schmidt-Thieme, 2010;Cong et al., 2015), the RESCAL model is a more suitable choice of tensor factorization for relational learning, because a more reasonable low-rank structure is used. More detail will be illustrated in the next subsection.

Related work in relational learning
Extensive methods have been explored to do link prediction with applications in all areas including knowledge graph completion Ji et al., 2016;Nguyen et al., 2017), biological network analysis (Almansoori et al., 2012;Sulaimany et al., 2018;Peng et al., 2020), and social network analysis (Al Hasan and Zaki, 2011;Liben-Nowell and Kleinberg, 2007;Almansoori et al., 2012). One line of the research is the translation-based models.
Among them, TransE (Bordes et al., 2013) is a very simple and effective one, where a relational fact is represented by a vector (h,r,t), which represents head, relationship, tail respectively. TransE assumes the score function f r (h, t) = h + r − t 2 2 , whose value is low if (h,r,t) holds and high otherwise. Later on, TransH  is proposed to overcome the flaws of TransE in dealing with 1-to-N, N-to-1 and N-to-N relation, where relations are modeled as a vector r on a hyperplane with a normal vector w r . The score r hw r , t ⊥ = t − w T r tw r and w r 2 = 1. Another algorithm TransR  follows the same idea except that they embed entities and relations in distinct entity space and relation space. A more complicated model TATEC (Garcia-Duran et al., 2016) encompasses previous embedding works by combining two-way interactions with high-capacity three-way ones, which achieves even higher accuracy. Other variants of the translation-based embedding models including TransM (Fan et al., 2014), TransA (Xiao et al., 2015), STransE (Nguyen et al., 2016) , and TransF (Feng et al., 2016). While these translation-based models are effective and showed good performance, only the observed triplets with a set of sampled unobserved triplets are taken into account during optimization.
Another important class of models represent relational facts as tensors, which draws a growing interest due to their natural representation of multi-relational data and its ability to utilize all the triplets' information in the data. Franz et al. (2009) proposed the TripleRank method, where they used PARAFAC tensor factorization to perform link prediction. The data is represented by a tensor X ∈ R k×ℓ×m , which is decomposed using X = n k=1 λ k · U k 1 • U k 2 • U k 3 , such that U 1 ∈ R k×n , U 2 ∈ R ℓ×n , U 3 ∈ R m×n . and "• ′′ denotes outer product. Bader et al. (2007) used a family of models called DEDICOM from the psychometrics literature to do tensor factorization in social networks, which shows the advantages of tensor factorization in dealing with such relational datasets. Nickel et al. (2011) proposed the RESCAL model that modified the DEDICOM by imposing a different constraint, which is more appropriate to use for multi-relational data.
Specifically the model is X k ≈ AW k A T where A contains the latent-component representation of the entities and W k models the interaction of the latent component in the k-th relationship. RESCAL in Nickel et al. (2011) uses the ASALSAN algorithm to find optimal parameters, where matrix A and matrices W k , k = 1, 2, ..., K are updated using efficient alternating least squares (ALS). The RESCAL model demonstrated its scalability by conducting experiments on YAGO 2 ontology data, which consists of 2.6 million entities, 33 million facts and 87 predicates. While the RESCAL achieves good performance in efficiency, scalability and accuracy, binary relational data are considered to be continuous in the model. As a result, an appropriate threshold is needed to be given before link prediction.
To address this issue, Nickel et al. (2013) proposed a logistic tensor factorization methods.
In this model, they assume each element in X k follows a Bernoulli distribution, which outperforms the original RESCAL model in terms of prediction accuracy measured by AUC.
But a serious drawback of the logistic tensor factorization methods is its poor scalability where dense matrices AW k A T needs to be computed in each iteration, whose complexity is quadratic in the number of entities. So it would be infeasible to use logistic tensor factorization in knowledge base like YAGO2. Some other extensions of RESCAL model that take the binary nature of the data into account have been explored by London et al. (2013) and Hu et al. (2016), Ghadermarzy et al. (2018), where the first one discussed different loss functions and used L-BFGS to estimate the optimal parameters. and the latter two are fully bayesian based factorization models. Our proposed PTF model is an improvement on both the RESCAL model and logistic tensor factorization models in that it inherits the computation efficiency from the former and enjoys the prediction accuracy from the latter.

Related work in network analysis
Our work is also closely related with researches of network analysis in statistics literature, especially the literatures where the data is represented as a multi-way tensor data. Two common multi-way tensor related topics have been studied extensively.
One important topic in the area is to describe the high-dimensional multi-way tensor data with some low-rank representations. A commonly adopted frameworks is a model of the form X = M + E, where X denotes the observed tensor data, E is noise and M is a mean array that is often assumed to be of low rank and can characterize the signal of interest with fewer parameters. To infer the low-dimensional representation of the connection patterns M, different types of models have been explored (Westveld et al., 2011;Hoff, 2011;Schein et al., 2015;Minhas et al., 2016;Hoff et al., 2016;Vervliet et al., 2019).
The other closely related line of research focuses on dynamic networks, which depict how the network data evolves over time. The goal is to predict future networks and capture the linkage evolution. It has been posited as a major topic of interest in many areas. A great number of researches within statistical frameworks have been proposed (Sarkar and Moore, 2005;Sarkar et al., 2007;Foulds et al., 2011;Xu and Hero, 2014;Durante and Dunson, 2014;Sewell and Chen, 2015;Hoff, 2011Hoff, , 2015, where different hidden structures are assumed. Despite advances in statistical methods of dynamic network analysis, deep models with the advantage in modelling non-linear transformation of network structure have also been explored (Perozzi et al., 2014;Wang et al., 2016;Grover and Leskovec, 2016;Trivedi et al., 2017). In this paper, while we are also modeling with three-way tensors, we specifically target at the task of missing link prediction and restrict to the third slice in the tensors to represent different types of relations between entities. We expect that for different types of relations, the low-rank representations of the entities will be interacted in different ways, which will be characterized by relation-specific latent matrices in our model.

Model Description
Tensors are denoted by bold upper case letters, such as X, Z and W. Denote the k-th slice of a tensor Z by Z k , the element in the i-th row, j-th column of Z k by z ijk . Matrices are denoted by capital letters such as A, whereas a i denotes i-th row of matrix A.

Model Framework
Assume that there are K types of relationships between N entities. The relational facts can be represented as a three-way tensor X ∈ {−1, 1} N ×N ×K with some missing values.
X k is an asymmetric N × N matrix representing the facts in the k-th type of relation.
Specifically, x ijk = 1 and x ijk = −1 stand for existing (valid) and non-existing (invalid) triplet (i-th entity, k-th relationship, j-th entity) respectively. The goal of the proposed model is to estimate P r(x ijk = 1) for all missing triplets given the observed relations.
A latent three-way tensor variable Z is introduced with the same dimension as X, where each element z ijk follows a Gaussian distribution. More specifically, for the k-th relation slice, Here r is the factorization rank, vec stands for the vectorization operator and I is the identity matrix.
The i-th row of A can be interpreted as the latent component representation for the i-th entity, and W k can be viewed as the interaction effects of latent components for the k-th relation. We expect that for different relations, the latent representations for entities will interact with each other in different ways, which can be captured by the different in W k .
Note that the bilinear product AW k A T can be written as A constraint that x ijk = sign(z ijk ) for all i, j, k is imposed, which implies the probit link so that P r(x ijk = 1) = P r(z ijk > 0). There are two reasons for using the probit link based on the latent tensor variable Z. First, the model outputs probabilities for each of the unknown triplets to be valid, which allows interpretable, probabilistic inference. Second, the use of the probit link enables an efficient iterative estimation method.

EM Algorithm
The basic idea of our model fitting is to maximize the log-likelihood, which is equivalent to minimizing cross-entropy loss. Since a latent variable is included in our model, directly estimating of model parameter would be numerically intensive. We use the classic EM algorithm to handle the latent variable.
EM algorithm is an iterative method to perform parameter estimation involving latent variables. It was first formally presented by Dempster et al. (1977) and has been a popular tool in numerical optimization. It consists of a E-step and a M-step.

E-step
The complete data log-likelihood, where we assume the latent variable Z is observed, can be written as where C is a constant. In addition, we impose the constraint such that sign(z ijk ) = x ijk for those observed x ijk , i.e. z ijk should be valid if x ijk = 1 and invalid if x ijk = −1.
Then, the Q function can be written as where C 2 is another constant such that and E A (t) ,W (t) denotes the expectation at the parameter values from last iteration t.

So in the M step, it is equivalent to minimize
Since the latent variable is assumed to follow a Gaussian distribution with constraint x ijk = sign(z ijk ), we have at iteration t such that where µ ijk = a T i W k a j , and a i , a j are i-th and j-th row from matrix A from last iteration. TN denotes truncated Gaussian distribution, with parameter space truncated to the domain defined after semicolon. Then, the expectation of the truncated Gaussian distribution can be computed as For simplicity of the notations, we denote the expectation matrix E A (t) ,W (t) (Z k |X k ) by E k . For k=1,...K, we define M k as sparse matrices such that where φ(.) and Φ(.) are the probability density function and cumulative distribution function for the standard Gaussian distribution.
Then we have E k = AW k A T + M k . At each iteration, we only need to compute the sparse matrices M k 's rather than the dense matrices E k 's. This technique will be illustrated in the M-step, which significantly saves the computation cost of the algorithm.

M-step
In the M-step, our goal is to maximize Q function, which is equivalent to solve the following optimization problem If E k 's were calculated in the E-step, then the update for A and W would be the same as RESCAL model. However, we only need the matrices M k 's to carry out the M-step.
Most of the elements in M k 's are usually zero, since most relational datasets are sparse, i.e. only a small fraction of relational facts are known. Not having to evaluate E k 's saves us O(N 2 r) computation cost in each iteration. The parameter updating rules are based on the alternating least-squares (ALS) algorithm as follows.

Update W in M Step:
When A is fixed, our objective function about W k is which becomes a least square problem.
Let β = A ⊗ A, where ⊗ is the Kronecker product.
Then the update rule is In the update for W k , the most expensive computation lies in computing (β T β) −1 β T .
To reduce the complexity, To reduce the complexity, we apply the singular value decomposition(SVD) to A as suggested by Chang et.al (2014), i.e. A = UΣV T with U and V being orthogonal matrices and Σ being a diagonal matrix. Then, The update rule can be rewritten as Update A in M Step: In order to simultaneously solve for the left matrix and right matrix A, we stack the data side by side as suggested by ASALSAN in Bader et al. (2007). LetĒ = (E 1 , E T 1 , ... , E K , E T K ) ,W = (W 1 , W T 1 , ... W K , W T K ), then it is equivalent to writeĒ = AW (I 2k ⊗ A T ). By setting its gradient to 0, we can get the update rule of A as We can compute E k AW T k and E T k AW k in the following way. Similarly,

Parameter Initialization
Note that the objective function is non-convex. The optimization process searches for locally optimal solutions and is not guaranteed to converge to a global optimal. Initialization is crucial for the performance. Similar to RESCAL (Nickel et al., 2011), the initialization process is performed with a singular value decomposition overX = k (X k +X T k ) such that X = UΣV T . The initial value of A is set to U. W is randomly initialized with each element independently drawn from a standard normal distribution. We find such a initialization work well during our simulation study and real data applications.

Complexity Analysis
To ensure the algorithm is tractable, it is important to study the computation complexity for the learning algorithm. Recall that N denotes the number of entities, K denotes the number of relationship types, r denotes the rank for the factorization. Let T denote the number of known triplets. The total number of triplets is N 2 K. When only a small fraction of the triplets are known to be valid / invalid, T is much smaller than N 2 K and the matrices M k 's in our algorithm are sparse matrices with only T non-zero elements.
The computation complexity is O(T r 2 ) for each E-step and O(NKr 2 + T Nr + Kr 3 ) for each M-step, which includes O(NKr 2 +T r+Kr 3 ) for updating A and O(Nr 2 +T Nr+r 3 ) for updating W. In total, for each EM iteration, we need O(NKr 2 +T Nr +Kr 3 ) computation, which is linear in N and K, T and cubic in r.
In terms of memory complexity, the tensors and matrices we need to keep in memory include a sparse three-way data X of size N × N × K with T non-zero elements, whose memory complexity is O(T ), an N × r dense matrix A, and a r × r × K dense tensor.
Besides, we also need to store a sparse tensor M that has the same dimension and sparsity as X.

Link Prediction
After the convergence criterion is achieved, we obtain an estimated matrixÂ and an estimated tensorŴ. Then the link prediction can be done in the following way.
For any triplet (i,j,k) that is unobserved in the data X, we have z ijk = a T i W k a j + ǫ by the definition of latent variable Z with ǫ ∼ N(0, 1), which implies i.e. the probability of the triplet (i,j,k) being valid is Φ(â T iŴ kâj ).

Simulation Setup
In this section, we demonstrate the performance of the proposed estimators using the simulated data. We compared the three factorization-based models, including the RESCAL model, the RESCAL Logit model and our PTF model. Simulation studies with different settings and two generating models are considered. Both the in-sample performance and the out-of-sample performance are evaluated. For the in-sample estimation, we use the canonical correlations between A andÂ to evaluate the parameter estimation. The reason to use the canonical correlations is that the bilinear product AW k A T can be written as The two probabilistic model, the RESCAL Logit model and our proposed PTF are used to generate the data. We first generate entries in A from N(0, 1). The elements in W k for k = 1, ...K are independently drawn from a Gaussian distribution with variance 1 and means µ k , which independently follow a Unif orm(−2, −1) distribution. The error term ǫ is generated from i.i.d N(0, 1).
• Model 1: The proposed Probit Tensor Factorization model, where for i = 1,2,...,N, j = 1,2,...,N, k = 1,... K • Model 2: The RESCAL Logit model  with Table1 summarizes our two simulation settings. For each simulated data, we randomly set 50%, 70% or 90% of links as missing, which are the triplets to make prediction for.  We compare our proposed model with two other tensor factorization models, which are the RESCAL model (Nickel et al., 2011) and RESCAL Logit model . The results are averaged over 100 repetitions. We set the maximum number of iterations at 500. The maximum iteration in the M step in our model is set to 3. For RESCAL, we choose a threshold for the output score such that the proportion of zeros and ones in testing match the proportions in training.

Simulation results
In Table 2, we summarized the performances of the three tensor factorization models when the data is generated from the PTF model. For both in-sample estimation of the latent component matrix A and the out-of-sample link prediction, the PTF model is the best among the three. While all of these models have good performance in in-sample estimation, the PTF model shows its advantage in out-of-sample link prediction. The results are similar in Table 3, where the RESCAL Logit is used as the generating model.

Correct in-sample estimation on the latent component matrices provides opportunities
for future quantitative analysis with the entities. For example, it can be used to cluster the entities. The out-of-sample AUC measures how well the models performed in link prediction, which is one of the most common tasks in relational learning. From the above two tables, we find that our PTF model achieve good results even when the data is generated from RESCAL Logit model. While the RESCAL Logit model achieved satisfactory prediction performances when the missing proportions are low, it loses its power as the missing    We used an information criteria proposed by Shi et al. (2019) for rank selection, which helps to select the correct rank most of the time during simulation. We are also interested in the model performance when the factorization rank is mis-specified in the PTF model.
The Table 5 summarized the average AUCs over 100 petitions when rank is mis-specified.
Here we set N to 200, K to 10 and missing proportion to 50%. Data are generated from the PTF model. Even when rank is mis-specified, the PTF model still achieves satisfactory performance.

Real data application
This section presents the results from experiments with three common small-and moderatesized benchmark datasets: UMLS, Kinships, Nations to evaluate the predictive performance of our model.

Datasets
In order to compare our model with current state-of-art relational learning models, the evaluation was carried out on the following three benchmark datasets: UMLS The data is gathered by McCray et al. (2001)

Experimental Results
We followed the experimental settings used by Kemp et al. (2006). For each of the three datasets, a ten-fold cross-validation was performed with each triplet as one statistical unit.
The triplets in the test set are treated as missing. The goal is to accurately predict the missing triplets. Since all the relational facts are known in the datasets, AUC is used to evaluate the prediction accuracy. During the experiment, the maximum number of iterations in the M step is 2 and the maximum iteration number is 500.
We consider different types of methods as our competitor. MRC Kok and Domingos (2007) is a Markov-logic based model. TransE (Bordes et al., 2013) and SME (Bordes et al., 2014) are two popular energy-based models. Results from tensor factorization related models are also reported, including LFM (Jenatton et al., 2012), RESCAL ALS (Nickel et al., 2011), RESCAL Logit , TATEC (Garcia-Duran et al., 2016), BPBFM2 (Hu et al., 2016). In addition, we include a statistical network analysis model Latentnet (Krivitsky and Handcock, 2008), which was applied to each slice of the tensor separately. Table 7 summarizes the results of the three experiments. We can see that our proposed and various order models perform much better than the two-way models. AUC from PTF, RESCAL Logit, RESCAL are all above 0.95 and the PTF performs the best. In Nations dataset, only LFM and PTF achieve AUC higher than 0.9 and the PTF model still performs the best. In summary, in terms of prediction accuracy in these three datasets, the proposed PTF performs uniformly better than the competing methods.

Discussion
In this paper, we proposed a probit tensor factorization model, which is an extension of tensor factorization model for binary data. The two parameters in our model, A and W, correspond to the low-rank representation of entities and characterization of different relations. We illustrate our model fitting procedure under the framework of EM algorithm.
Scalability issue is addressed by showing that the complexity of our algorithm is linear in the number of entities, the number of known facts and the number of relationship types.
We conducted experiments on three common medium sized datasets, where our method obtains uniformly better results than the other state-of-art models. One major advantage of our model is that it adopts the open world assumption (OWA) and allows a natural use of invalid facts in model. And it outputs an interpretable probabilistic result that can be further used in quantitative analysis.
While our model is designed to deal with asymmetric relations, it can also be used to perform link predictions for symmetric relations. We need to impose a constraint that the latent representation for relations W k 's are symmetric metric. With the constraint, the update rule for W k can be formulated as a quadratic programming problem.
Another interesting topic is how to decide the factorization rank, where both computation cost and predictive performance should be taken into consideration. In the literature, Mazumder et al. (2010) applied regularization to select the rank for singular value decomposition, Nickel et al. (2011) suggested to use cross-validation for rank selection in RESCAL model and Sun et al. (2017) proposed to use Bayesian information criteria (BIC) (Schwarz et al., 1978) for sparse CP decomposition. Recently, Shi et al. (2019) proposed a general class of information criteria to determine the rank for RESCAL model and showed their method achieved good finite sample properties. The proposed method is also applicable to our model and we found it works well in our simulation.
The PTF model can be extended in several directions. One way is to incorporate prior information. For example, by utilizing prior domain knowledge of type constraints that is suggested by Chang et al. (2014), we can further improve the efficiency. And more prior information, e.g. logical structure of the relations, attributes of the entities and so on are expected to be included in tensor factorization model to further enhance the model performance.

A Appendix
A.1 Runtime for the EM algorithm We summarized runtime for both E step and M step with different entity size N, relation size K, number of observed triplets T and decomposition rank r. Note that the computation complexity is O(T r 2 ) for each E-step and O(NKr 2 +T r 2 ) for each M-step. All runtimes are averaged over 100 repetitions. The runtime performance is aligned with our expectation.