Relational graph convolutional networks: a closer look

Message passing

We will begin by describing the basic Graph Convolutional Network (GCN) layer for directed graphs¹ . This will serve as the basis for the RGCN in “Extending GCNs for multiple relations”.

The GCN (Kipf & Welling, 2017) is a graph-to-graph layer that takes a set of vectors representing nodes as input, together with the structure of the graph and generates a new collection of representations for nodes in the graph. A directed graph is defined as $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is a set of vertices (nodes) and $⟨i,j⟩\in \mathcal{E}$ is a set of tuples indicating the presence of directed edges, pointing from node i to j. Equation (1) shows the message passing rule of a single layer GCN for an undirected graph, $\mathcal{G}$.

(1) $$H=\sigma \left(AXW\right),$$

Here, X is a node feature matrix, W represents the weight parameters, and $\sigma$ is a non-linear activation function. A is a matrix computed by row-normalizing² the adjacency matrix of the graph $\mathcal{G}$. The row normalization ensures that the scale of the node feature vectors do not change significantly during message passing. The node feature matrix, X, indicate the presence or absence of a particular feature on a node.

Typically, more than a single convolutional layer is required to capture the complexity of large graphs. In these cases, the RGCN layers are stacked one after another so that the output of the preceeding RGCN layer $H^{(l-1)}$ is used as the input for the current layer $H^{(l)}$, as shown in Eq. (2). In our work, we will use superscript l to denote the current layer.

(2) $$H^{(l)}=\sigma \left(A{H}^{(l-1)}W\right),$$

If the data comes with a feature vector for each node, these can be used as the input X for the first layer of the model. If feature vectors are not available, one-hot vectors, of length N with the non-zero element indicating the node index, are often used. In this case, the input X becomes the identity matrix I, which can then be removed from Eq. (1).

We can rewrite Eq. (1) to make it explicit how the node representations are updated based on a node’s neighbors:

(3) $${\mathbf{h}}_i=\sigma \left[\sum _{j\in N(i)}{\displaystyle \frac{1}{|N(i)|}}{{\mathbf{x}}_i}^{T}W\right].$$

Here, ${\mathbf{x}}_i$ is an input vector representing node i, ${\mathbf{h}}_i$ is the output vector for node i, and N(i) is the collection of the incoming neighbors of i, that is the nodes j for which there is an edge $⟨$j, i $⟩$ in the graph. For simplicity, the bias term is left out of the notation but it is usually included. We see that the GCN takes the average of i’s neighbouring nodes, and then applies a weight matrix W and an activation $\sigma$ to the result. Multiplying ${{\mathbf{x}}_i}^{T}W$ by ${\displaystyle \frac{1}{|N_i|}}$ means that we sum up all the feature vectors of all neighboring nodes. This makes every convolution layer permutation equivariant, and that is: if either the nodes (in A) are permuted, then the output representations are permuted in the same way. Overall, this operation has the effect of passing information about neighboring nodes to the node of interest, i and is called message passing. Message passing is graphically represented in Fig. 1 for an undirected graph, where messages from neighboring nodes (a − e) are combined to generate a representation for node i. After message passing, the new representation of node i is a mixture of the vector embeddings of neighboring nodes.

If a graph is sparsely connected, a single graph-convolution layer may suffice for a given downstream task. Using more convolutional layers encourages mixing with nodes more than 1-hop away, however it can also lead to output features being oversmoothed (Li, Han & Wu, 2018). This is an issue as the embeddings for different nodes may be indistinguishable from each other, which is not desirable.

In summary, GCNs perform the following two operations: (1) They replace each node representation by the average of its neighbors, and (2) they apply a linear layer with a nonlinear activation function $\sigma$. There are two issues with this definition of the GCN. First, the input representation of node i does not affect the output representation, unless the graph contains a self-loop for i. This is often solved by adding self-loops explicitly to all nodes. Second, only the representations of nodes that have incoming links to i are used in the new representation of i. In the relational setting, we can solve both problems elegantly by adding relations to the graph which we will describe in the next section.

Extending GCNs for multiple relations

In this section, we explain how the basic message passing framework can be extended to relational graphs, such as Knowledge Graphs. We define a Knowledge Graph as a directed graph with labelled vertices and edges. Formally, this graph can be defined as $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{R})$, where $\mathcal{R}$ represents the set of edge labels (relations) and $⟨s,r,o⟩\in \mathcal{E}$ is a set of tuples representing that a subject entity s and an object entity o are connected by the relation $r\in \mathcal{R}$.

The Relational Graph Convolutional Network extends graph convolutions to Knowledge Graphs by accounting for the directions of the edges and handling message passing for different relations separately. Equation (4) is an extension of the regular message passing rule (Eq. (1)).

(4) $$H=\sigma \left(\sum _{r=1}^R{A}_{r}X{W}_r\right),$$where R is the number of relations, $A_r$ is an adjacency matrix describing the edge connection for a given relation r and $W_r$ is a relation-specific weight matrix. The extended message passing rule defines how the information should be mixed together with neighboring nodes in a relational graph. In the message passing step, the embedding is summed over the different relations³ . We can rewrite Eq. (4) to show how the node representations are updated based a node’s neighbors connected via different relations:

(5) $${\mathbf{h}}_i=\sigma \left[\sum _r^R\sum _{j\in N_r(i)}{\displaystyle \frac{1}{|N_r(i)|}}{\mathbf{x}}_i^T{W}_r\right],$$where $N_r(i)$ is the collection of the incoming neighbors of i with the relation r.

With the message passing rule discussed thus far, the problem is that for a given triple $⟨$s, r, o $⟩$ a message is passed from s to o, but not from o to s. For instance, for the triple $⟨$Amsterdam, located_in, The_Netherlands $⟩$ it would be desirable to update both Amsterdam with information from The_Netherlands, and The_Netherlands with information from Amsterdam, while modelling the two directions as meaning different things. To allow the model to pass messages in two directions, the graph is amended inside the RGCN layer by including inverse edges: for each existing edge $⟨$s, r, o $⟩$, a new edge $⟨o,r^{\mathrm{\prime }},s⟩$ is added where $r^{\mathrm{\prime }}$ is a new relation representing the inverse of r. A second problem with the naive implementation of the (R)GCN is that the output representation for a node i does not retain any of the information from the input representation. To allow such information to be retained, a self-loop $⟨s,r_s,s⟩$ is added to each node, where $r_s$ is a new relation that expresses identity. Altogether, if the input graph contains R relations, the amended graph contains 2R + 1 relations: ${\mathcal{R}}^{+}=\mathcal{R}\cup {\mathcal{R}}^{\mathrm{\prime }}\cup {\mathcal{R}}_s$. This is graphically represented in Fig. 2 for a directed graph.

Reducing the number of parameters

We use $N_{in}$ and $N_{out}$ to represent the input and output dimensions of a layer, respectively. While the GCN (Kipf & Welling, 2017) requires $N_{in}\times N_{out}$ parameters, relational message passing uses $R^{+}\times N_{in}\times N_{out}$ parameters. In addition to the extra parameters required for a separate GCN for every relation, we also face the problem that Knowledge Graphs do not usually come with a feature vector representing each node. As a result, as we saw in the previous section, the first layer of an RGCN model is often fed with a one-hot vector for each node. This means that for the first layer $N_{in}$ is equal to the number of nodes in the graph.

In their work, Schlichtkrull et al. (2018) introduced two weight regularisation techniques: (1) Basis Decomposition and (2) Block Diagonal Decomposition. We believe that these matrix decomposition techniques help to reduce the number of parameters (i.e., increasing parameter efficiency) and thus, prevent the model from overfitting to the training datasets. Figure 3 shows visually how the two different regularisation techniques work.

Basis Decomposition does not create a separate weight matrix $W_r$ for every relation. Instead, the matrices $W_r$ are derived as linear combinations of a smaller set of B basis matrices $V_b$, which is shared across all relations. Each matrix $W_r$ is then a weighted sum of the basis vectors with component weight $C_{rb}$:

(6) $$W_r=\sum _{b=1}^B{C}_{rb}{V}_b,$$

Both the component weights and the basis matrices are learnable parameters of the model, and in total they contain fewer parameters than $W_r$. With lower number of basis functions, B, the model will have reduced degrees of freedom and possibly better generalisation.

Block Diagonal Decomposition creates a weight matrix for each relation, $W_r$, by partitioning $W_r$ into ${\displaystyle \frac{N_{in}}{B}}$ by ${\displaystyle \frac{N_{out}}{B}}$ blocks, and then fixing the off-diagonal blocks to zeros (shown in Fig. 3)⁴ . This deactivates the off-diagonal blocks, such that only the diagonal blocks are updated during training. An important requirement for this decomposition method is that the width/height of $W_r$ need to be divisible by B.

(7) $$W_r=\left[\begin{array}{cccc}{Q}_{r1}& 0& \dots & 0\\ 0& Q_{r2}& 0& 0\\ ⋮& 0& \ddots & ⋮\\ 0& 0& \dots & Q_{rB}\end{array}\right]$$

Here, B represents the number of blocks that $W_r$ is decomposed into and $Q_{rb}$ are the diagonal blocks containing the relation-specific weight parameters. Equation (7) shows that taking the direct sum of $Q_r$ over all the blocks gives $W_r$, which can also be expressed as the sum of $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(Q_{rb})$ over all the diagonal elements in the block matrix. The higher the number of blocks b, the lower the number of trainable weight parameters for each relation, $W_r$, and vice versa. Block diagonal decomposition is not applied to the weight matrix of the identity relation $r_s$, which we introduced in “Extending GCNs for multiple relations” to add self-loops to the graph.

Einstein summation

Message passing requires manipulating high-dimensional matrices and tensors using many different operations (e.g., transposing, summing, matrix-matrix multiplication, tensor contraction). We will use Einstein summation to express these operations concisely.

Einstein summation (einsum) is a notational convention for simplifying tensor operations (Kuptsov, 2022). Einsum takes two arguments: (1) an equation⁶ in which characters identifying tensor dimensions provide instructions on how the tensors should be transformed and reshaped, and (2) a set of tensors that need to be transformed or reshaped. For example, einsum ( $ik,jk\to ij,A,B$) represents the following matrix operation:

(8) $$C_{ij}=\sum _k{A}_{ik}\cdot B_{jk},$$

The general rules of an einsum operation are that indices which are excluded from the result are summed out, and indices which are included in all terms are treated as the batch dimension. We use einsum operations in our implementation to simplify the message passing operations.

Sparsity

Since many graphs are sparsely connected, their adjacency matrices can be efficiently stored on memory as sparse tensors. Sparse tensors are memory efficient because, unlike dense tensors, sparse tensors only store non-zero values and their indices. We make use of sparse matrix multiplications⁷ . For sparse matrix operations on GPUs, the only multiplication operation that is commonly available is multiplication of a sparse matrix S by a dense matrix D, resulting in a dense matrix. We will refer to this operation as spmm(S, D). For our implementation, we endeavour to express the sparse part of the RGCN message passing operation (Eq. (4)), including the sum over relations, in a single sparse matrix multiplication.

Stacking trick

Using nested loops to iteratively pass messages between all neighboring nodes in a large graph would be very inefficient. Instead, we use a trick to efficiently perform message passing for all relations in parallel.

Edge connectivity in a relational graph is represented as a three-dimensional adjacency tensor $A\in {\mathbb{R}}^{R^{+}\times N\times N}$, where N represents the number of nodes and $R^{+}$ represents the number of relations. Typically, message passing is performed using batch matrix multiplications as shown in Eq. (4). However, at the time of writing, batch matrix operations for sparse tensors are not available in most Deep Learning libraries. Using spmm is the only efficient operation available, so we stack adjacency matrices and implement the whole RGCN in terms of this operation.

We augment A by stacking the adjacency matrices corresponding to the different relations $A_r$ vertically and horizontally into $A_v\in {\mathbb{R}}^{(N+R^{+})\times N}$ and $A_h\in {\mathbb{R}}^{N\times (N+R^{+})}$, respectively.

(9) $$A_v=\left[\begin{array}{c}{A}_1\\ A_2\\ ⋮\\ A_{1^{\mathrm{\prime }}}\\ A_{2^{\mathrm{\prime }}}\\ ⋮\\ A_s\end{array}\right]$$

(10) $$A_h=\left[\begin{array}{c}{A}_1{A}_2\cdots A_{1^{\mathrm{\prime }}}{A}_{2^{\mathrm{\prime }}}\cdots A_s\end{array}\right]$$

Here, $[\cdot ]$ represents a concatenation operation, and $A_v$ and $A_h$ are both sparse matrices. By stacking $A_r$ either horizontally or vertically, we can perform message passing using sparse matrix multiplications rather than expensive dense tensor multiplications. Thus, this trick helps to keep the memory usage low.

Algorithm 1 shows how message passing is performed using a series of matrix operations. All these are implementations of the same operation, but with different complexities depending on the shape of the input. In Theorem 1 (see the Appendix), we prove that the the stacked formulation is equivalent to the original RGCN formulation.

1. If the inputs X to the RGCN layer are one-hot vectors, X can be removed from the multiplication. The featureless message passing simply multiplies A with W, because the node feature matrix X is not given. Note that X, in this case, can also be modelled using an identity matrix I. However, since AW = AIW, we skip this step to reduce computational overhead.

2. In the horizontal stacking approach, X multiplied with W. This yields the XW tensor, which is then reshaped into a $N{R}^{+}\times N$ matrix. The reshaped XW matrix is then multiplied with $A_h$ using spmm.

3. In the vertical stacking approach, the X is mixed with $A_v$ using spmm. The product is reshaped into a tensor of dimension $R^{+}\times N\times N$. The tensor AX is then multiplied with W.

Any dense/dense tensor operations is implemented with the einsum operation. For reasons highlighted in “Sparsity” it is desirable to construct $A_r$ as a sparse tensor. However, in PyTorch (Paszke et al., 2019) sparse/dense operations only allow multiplication of sparse matrix by dense matrix. We accept this as a limitation of the current version of PyTorch. To circumvent this issue, we use the stacking trick to construct $A_r$ as a sparse matrix and use the spmm (S, D) operation. Thus, the stacking trick is a memory efficient operation for the message passing operation in the RGCN layers.

Any dense/dense tensor operations can be implemented with einsum, but in the current version of PyTorch sparse/dense operations only allow multiplication of sparse matrix by dense matrix. Therefore, the stacking trick provides a memory efficient way to perform message passing operation in the RGCN layers (see Sparsity and Stacking trick). The vertical stacking approach is suitable for low dimensional input and high dimensional output, because the projection to low dimensions is done first. While the horizontal stacking approach is good for high dimensional input and low dimensional output as the projection to high dimension is done last. These matrix operations are visually illustrated in Fig. 4.

Thus far, we focused on how Relational Graph Convolutional layers work and how to implement them. As mentioned in “Literature review”, RGCNs can be used for many downstream tasks. Now, we will discuss how these graph convolutional layers can be used as building blocks in larger neural networks to solve two downstream tasks implemented in the original RGCN article (Schlichtkrull et al., 2018): node classification and link prediction. In the next two sections, we detail the model setup, our reproduction experiments and new configurations of the models. We begin with node classification.

Model setup

Figure 5 shows the node classification models described in (Schlichtkrull et al., 2018), using a two-layer model. Full-batch training is used for training the node classification model, meaning that the whole graph is represented in the input adjacency matrix A for the RGCN. The input is the unlabeled graph, the output are the class predictions and the true predictions are used to train the model. The first layer of the RGCN is ReLU activated and it embeds the relational graph to produce low-dimensional node embeddings. The second RGCN layer further mixes the node embeddings. Using softmax activation the second layer generates a matrix with the class probabilities, $Y\in {\mathbb{R}}^{N\times C}$, and the most probable classes are selected for each unlabelled node in the graph. The model is trained by optimizing the categorical cross entropy loss:

(11) $$\mathcal{L}=-\sum _{i=1}^{\mathcal{Y}}\sum _{c=1}^C{t}_{ic}\ln {\mathbf{h}}_{ic}^{(L)},$$where $\mathcal{Y}$ is the set of labelled nodes, C represents the number of classes, $t_{ic}$ is one-hot encoded ground truth labels and ${\mathbf{h}}_{ic}^{(L)}$ represents node representations from the RGCN. The last layer (L) of the RGCN is softmax-activated. The trained model can be used to infer the classes of unlabelled nodes.

e-RGCN

In GCNs (Kipf & Welling, 2017), the node features X are represented by a matrix $X\in {\mathbb{R}}^{N\times F}$, where N is the number of nodes and F is the number of node features. When node features are not available, one-hot vectors can be used instead. An alternative approach would be to represent the features with continuous values $E\in {\mathbb{R}}^{N\times D}$, where D is the node embedding dimension.

In the GCN setting (Kipf & Welling, 2017), using one hot vectors is functionally very similar to using embedding vectors: the multiplication of the one hot vector by the first weight matrix W, essentially selects a row of W, which then functions as an embedding of that node. In the RGCN setting, the same holds, but we have a separate weight matrix for each relation, so using one-hot vectors is similar to defining a separate node embedding for each relation. When we feed the RGCN a single node embedding for each node instead, we should increase the embedding dimension D to compensate.

Initial experiments showed that simply replacing the node features X with node embeddings E⁸ results in a drop in performance on the node classification benchmark data. We believe that the additional parameters from E gave the model more flexibility to overfit to training data and thus leading to an overall performance drop on the test dataset. After some experimentation, we ended up with the following model, which we call the embedding-RGCN (e-RGCN). Its message passing rule is described in Eq. (12). The weight matrix is restricted to a diagonal matrix (with all off diagonal elements fixed to zero)⁹ and then the product is multiplied by the adjacency matrix.

(12) $$h=\sigma \left(\sum _{r=1}^R{A}_{r}E\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{w}}_r)\right),$$

Here, E is the node embeddings broadcasted across all the relations $\mathcal{R}$, $w_r$ is a vector containing weight parameters for relation r. Here, diag(⋅) is a function that takes a vector $\mathbf{l}\in {\mathbb{R}}^Q$ as an input and outputs a diagonal matrix $N\in {\mathbb{R}}^{Q\times Q}$, where the diagonal elements are elements from the original vector L.

Using a diagonal weight matrix improves parameter efficiency, while enabling distinction between relations. We created a new node classification model, where the first layer is an e-RGCN layer and the second layer is a standard RGCN (without regularisation) that predicts class probabilities. This model provides competitive performance with the RGCN, using only 8% of the parameters.

Node classification experiments

Model setup

We follow the procedure outlined by Schlichtkrull et al. (2018). Figure 6 shows a schematic representation of the link prediction model as described in the original article. During training, traditional link predictors (Bordes et al., 2013; Yang et al., 2015) simultaneously update node representations and learn a scoring function (decoder) that predicts the likelihood of the correct triple being true. RGCN-based link predictors introduce additional steps upstream.

We begin by sampling 30,000 edges from the graph using an approach called neighborhood edge sampling (see Edge sampling). Then, for each triple we generate 10 negative training examples, generating a batch size of 330,000 edges in total. Node embeddings $E\in {\mathbb{R}}^{N\times D}$ are generated using a standard embedding layer and are used as an input for the RGCN¹³ . The RGCN performs message passing over the sampled edges and generates mixed node embeddings. Finally, the DistMult scoring function (Yang et al., 2015) uses the mixed node embeddings and relation embeddings to compute the likelihood of a link existing between a pair of nodes. For a given triple $⟨$s, r, o $⟩$, the model is trained by scoring all potential combination of the triple using the function:

(13) $$f(s,r,o)=\sum _i{\left({\mathbf{e}}_s\mathbf{r}{\mathbf{e}}_o\right)}_i=x$$

Here, ${\mathbf{e}}_s$ and ${\mathbf{e}}_o$ are the corresponding node embedding of entities s and o, generated by the RGCN encoder. r is a low-dimensional vector of relation r, which is generated by the DistMult decoder. Relation embeddings $R\in {\mathbb{R}}^{R^{+}\times D}$ are generated using a standard embedding layer. As Schlichtkrull et al. (2018) highlighted in their work, the DistMult decoder can be replaced by any Knowledge Graph scoring function.

Similar to previous work on link prediction (Yang et al., 2015), the model is trained using negative training examples. For each observed (positive) triple in the training set, we randomly sample 10 negative triples (i.e., we use a negative sampling rate of 10). These samples are produced by randomly corrupting either the subject or the object of the positive example (the probability of corrupting the subject is 50%). Binary cross entropy loss¹⁴ is used as the optimization objective to push the model to score observable triples higher than the negative ones:

(14) $$\begin{array}{c}\hfill \mathcal{L}=-\sum _{(s,r,o,y)\in \mathcal{T}}y\log l(f(s,r,o))+(1-y)\log (1-l(f(s,r,o))),\end{array}$$where $\mathcal{T}$ is the total set of positive and negative triples, l is the logistic sigmoid function, and y is an indicator set to y = 1 for positive triples and y = 0 for negative triples. f(s, r, o) includes entity embeddings from the RGCN encoder and relations embeddings from the DistMult decoder.

Link prediction experiments

As in the original article, the models are evaluated using Mean Reciprocal Rank (MRR) and Hits@k (k = 1, 3 or 10). The Torch-RGCN model was trained for 7,000 epochs without early stopping¹⁶ . We monitored the training by evaluating the model at regular intervals (every 500 epochs) using the validation dataset. Schlichtkrull initialisation (see Appendix) was used to initialise all parameters in the link prediction models and in our reproductions. Schlichtkrull et al. (2018) trained their models on the CPU. Our Torch-RGCN implementation can be trained using GPU acceleration. Early stopping was not used. We used the hyperparameters described in (Schlichtkrull et al., 2018).

c-RGCN

The link prediction architecture presented in Schlichtkrull et al. (2018) does not represent a realistic competitor for the state of the art and is very costly to use on large graphs. Furthermore, a problem with the original RGCN link predictor is that we need high dimensional node representations to be competitive with traditional link predictors, such as DistMult (Yang et al., 2015), but the RGCN is expensive for high dimensions. However, we do believe that the basic idea of message passing is worth exploring further in a link prediction setting.

To show that there is promise in this direction, we offer a simplified link prediction architecture that uses a fraction of the parameters of the original implementation (Schlichtkrull et al., 2018) uses. This variant places a bottleneck architecture around the RGCN in the link prediction model, such that the node embedding matrix E is projected down to a lower dimension, C, and then the RGCN performs message propagation using the compressed node embeddings. Finally, the output is projected up back to the original dimension, D, and computes the DistMult score from the resulting high-dimensional node representations. We call this encoding network the compression-RGCN (c-RGCN). Figure 7 shows a schematic representation of the c-RGCN model. Equations (15) and (16) show the message passing rule for the first and second layer of the c-RGCN encoder, respectively. We selected a node embedding size of 128 and compressed it to a vector dimension of 16. We also include a residual connection by including E in the second layer of the c-RGCN. The residual connection allows the model, in principle, to revert back to DistMult if the message passing adds no value. If all RGCN weights are set to 0, we recover the original DistMult. The c-RGCN model can be trained on the GPU.

(15) $$H^1=\sigma \left(\sum _{r=1}^{R^{+}}A{W}_r{f}_{\theta }(E)\right),$$where $f(X)=X{W}_{\theta }+b$ with $W_{\theta }\in {\mathbb{R}}^{C\times D}$.

(16) $$H^2=E+g_{\varphi }\left(\sigma \left(\sum _{r=1}^{R^{+}}A{W}_r{H}^1\right)\right),$$where $g(X)=X{W}_{\varphi }+b$ with $W_{\varphi }\in {\mathbb{R}}^{D\times C}$.

As shown in Table 5, the c-RGCN does not perform much worse than the original implementation. However, it is much faster, and memory efficient enough for full batch evaluation on the GPU. There is a clear trade-off between compression size of the node embeddings and the performance in link prediction. While this result is far from a state of the art model, it serves as a proof-of-concept that there may be ways to configure RGCN models for a better performance/efficiency tradeoff.

Implications for RGCN usage

We believe that Relational Graph Convolutional Networks are still very relevant because it is one of the simplest members of the message passing models and is a good starting place for exploration of machine learning for Knowledge Graphs.

RGCNs clearly perform well on node classification tasks because the task of classifying nodes benefits from message passing. This means that a class for a particular node is selected by reasoning about the classes of neighboring nodes. For example, a researcher can be categorised into a research domain by reasoning about information regarding their research group and close collaborators. There is no direct way to perform node classification using traditional Knowledge Graph Embeddings (KGE) models, such as TransE and DistMult.

While the RGCN is a promising framework, in its current setting we found that the link prediction model proposed by Schlichtkrull et al. (2018) is not competitive with current state of the art (Ruffinelli, Broscheit & Gemulla, 2020) and the model is too expensive with considerably lower performance. In our article, we clarify that RGCN-based link predictors are extensions of KGE models (Ruffinelli, Broscheit & Gemulla, 2020), thus training RGCN to predict links will always be more expensive than using a state of the art KGE model. RGCN-based link predictor take several days to train, while state of the art relation models run in well under an hour (Ruffinelli, Broscheit & Gemulla, 2020).

To aid the usage of RGCN, we presented two new configurations of the RGCN:

e-RGCN. We propose a new variant of the node classification model which uses significantly less parameters by exploiting a diagonal weight matrix. Our results indicate that e-RGCN has the potential to perform competitively with the original model and deserves further investigation. The potential advantage of e-RGCN is that it can operate on graphs that are much larger than the benchmark datasets we used in our reproduction studies.

c-RGCN. We also present a proof-of-concept model that performs message passing over compressed graph inputs and thus, improves the parameter efficiency for link prediction. The c-RGCN has several advantages over the regular RGCN link predictor: (1) c-RGCN does not require sampling edges, because it is able to process the entire graph in full-batch, (2) c-RGCN takes a fraction of the time it takes to train an RGCN link predictor, (3) c-RGCN uses fewer parameters, and (4) it is straightforward to implement. Although the results for the c-RGCN are not as strong, this sets a path for further development towards efficient message models for relational graphs.

Reproduction

Evolving technologies pose several challenges for the reproducibilty of research artifacts. This includes frequent updates being made to existing frameworks, such as PyTorch and TensorFlow, often breaking backward compatibility. We were in a strong position to execute this reproduction: (1) an author of this article also worked on the original article, (2) we contacted one of the lead authors of this article who was very responsive and (3) we were able to run the original source code (https://github.com/tkipf/relational-gcn and https://github.com/MichSchli/RelationPrediction) inside a virtual environment. Nevertheless, we found it considerably challenging to make a complete reproduction. To explain why and to contribute to avoiding such situations in the future, we briefly outline the lessons we have learned during the reproduction.

Parameter Statistics. There were discrepancies between the description of the link prediction model in Schlichtkrull et al. (2018) and the source code. The source code reproduces the values similar to the MRR scores reported in Schlichtkrull et al. (2018). Thus, to reproduce the results we had to perform a deep investigation of the source code. Using the original source, we relied on comparing the parameter statistics and tensor sizes at various points in both models. Since these statistics are helpful to verify the correctness of an implementation, we believe this is a useful practice in aiding reproduction. For complex models with long runtimes, an overview of descriptive statistics of parameter and output tensors for the first forward pass can help to check implementation without running full experiments. We are publishing statistics for intermediate products that we obtained for the link prediction models (see Table 4).

Small dataset. We found that the link prediction datasets used by Schlichtkrull et al. (2018) were large and thus, impractical for debugging RGCN because it is costly to train them on large graphs. Using a smaller dataset (FB-Toy (Ruffinelli, Broscheit & Gemulla, 2020)) would enable quicker testing with less memory consumption. Thus, we report link prediction results on the FB-Toy dataset (see Table A2).

Training times. The training times were variable and strongly depended on the size of the graph, the number of relations and the number of epochs. Schlichtkrull et al. (2018) reported the computational complexity, but not practical training times. It turns out that this is an important source of uncertainty in verifying whether re-implementations are correct. We measured the runtimes, which includes training the model and using the pre-trained model for making inference. For 7,000 epochs, the link prediction runtimes for Torch-RGCN and c-RGCN on the WN18 dataset are 2,407 and 53 min, respectively. Node classification experiments took a few minutes to complete, because they only required 50–100 epochs. We encourage authors to report such concrete training times.

Hyperparameter Search. We found that hyperparameters reflect the complexity of the individual datasets. For example, AIFB, the smallest dataset, was not prone to overfitting. Whereas, the larger AM dateset required basis decomposition and needs a reduced hidden layer size. For link prediction, we were unable to identify the optimum hyperparameters for WN18, FB15k and FB15k-237 due to the sheer size of the hyperparameter space and long training times. We provide a detailed list of hyperparameter we use in our reproduction. While this is becoming more common in the literature, this serves as further evidence of the importance of this detailed hyperparameter reporting.

Other factors. We still faced the common challenges in software reproduction that others have long noted (Fokkens et al., 2013), including missing dependencies, outdated source code, and changing libraries. An additional challenge with machine learning models is that hardware (e.g., GPUs) now also can impact the performance of the model itself. For instance, while we were able to run the original link prediction code in TensorFlow 1.4, the models no longer seemed to benefit from the available modern GPUs. Authors should be mindful that even if legacy code remains executable for a long time, executing it efficiently on modern hardware may stop being possible much sooner. Here too, reporting results on small-scale experiments can help to test reproductions without the benefit of hardware acceleration.

Proofs

Definition 1. We define a relational graph $\mathcal{G}$ as a directed graph with labelled nodes and edges. F is the number of node features, R is the number of relations in the graph and N is the number of nodes in the graph. $A^r\in {\mathbb{R}}^{R\times N\times N}$ is a tensor that describes the edge connectivity for every relation. $A^h=\left[\begin{array}{c}{A}_1\\ A_2\\ ⋮\\ A_R\end{array}\right]$ is a matrix where the adjacency matrices corresponding to different relations are vertically stacked, and $A^v=\left[A_1{A}_2\cdots A_R\right]$ is a matrix where the adjacency matrices corresponding to different relations are horizontally stacked. $X\in {\mathbb{R}}^{N\times F}$ is a feature matrix. $W^r$ is a relation-specific weight matrix.

Theorem 1. The vertically and horizontally stacked versions of the RGCN are equivalent to the original definition.

Proof. To show that the stacking trick yields the same operation as the original RGCN definition.

First, we note that the non-linearity $\sigma$ is an element-wise operation. We can ignore this, and show that the input to the non-linearity is the same in all three cases.

To show this, we first express all three definitions in terms of the individual matrix elements. Writing the original R-GCN as

$G=\sum _{r=1}^R{A}^r\cdot X\cdot W^r$where X is the matrix of inputs, we can rewrite this operation in terms of the individual elements of G as

$$\begin{array}{c}{G}_{r,i,j}^{\mathrm{\prime }}=\sum _{k=1}^N{A}_{r,i,k}^r\cdot X_{r,k,j}\end{array}$$

$G_{r,i,j}^{\mathrm{\prime }\mathrm{\prime }}\begin{array}{c}=\sum _{k=1}^F{G^{\mathrm{\prime }}}_{r,i,k}\cdot W_{r,k,j}\end{array}$

$G_{r,i,j}\begin{array}{c}=\end{array}\sum _{k=1}^F{G}_{r,i,j}^{\mathrm{\prime }\mathrm{\prime }}$

Here, the first line shows a basic matrix multiplication, expressed element-wise $((AB{)}_{i,j}=\sum _k{A}_{i,k}\cdot B_{k,j}$, the second shows another matrix multiplication, and the third sums the result over all relations.

For the vertically stacked RGCN, we can write the computation of the output matrix V as

$\begin{array}{cc}{V^{\mathrm{\prime }}}_{i,j}& =\sum _{k=1}^N{A}_{i,k}^v\cdot X_{k,j}\\ {V^{\mathrm{\prime }\mathrm{\prime }}}_{r,i,j}& ={V^{\mathrm{\prime }}}_{(r-1)n+i,j}\\ {V^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}}_{r,i,j}& =\sum _{k=1}^F{V^{\mathrm{\prime }\mathrm{\prime }}}_{r,i,j}\cdot W_{r,k,j}\\ V_{i,j}& =\sum _{r=1}^R{V^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}}_{r,i,j}.\end{array}$

The first line again shows a matrix multiplication, but this time of the vertically stacked adjacency matrix with the feature matrix X. The second line reshapes this into a three-tensor separating our the relations. The third line is again a standard matrix multiplication and the final line sums out the relation.

We note that $V_{r,i,j}^{\mathrm{\prime }\mathrm{\prime }}=\sum _{k=1}^F{A}_{(r-1)n+i,k}^v{X}_{k,j}=\sum _{k=1}^F{A}_{r,i,k}^r{X}_{k,j}=G_{i,j}^{\mathrm{\prime }}$. It follows that V‴ = G″and V = G.

For the horizontally stacked RGCN, we can write the computation of the output matrix H as

$$\begin{array}{ll}{H^{\mathrm{\prime }}}_{r,i,j}& =\sum _{k=1}^F{X}_{r,i,k}\cdot W_{r,k,j}\\ {H^{\mathrm{\prime }\mathrm{\prime }}}_{(r-1)n+i,j}& ={H^{\mathrm{\prime }}}_{r,i,j}\\ H_{i,j}& =\sum _{k=1}^{RN}{A}_{i,k}^h\cdot {H^{\mathrm{\prime }\mathrm{\prime }}}_{k,j}.\end{array}$$

The first line shows a matrix multiplication of the input features by the weight matrix, per relation. The second line shows the reshaping of the result into a vertical stack of matrices and the final line shows a single matrix multiplication between the horizontally stacked adjacency matrix and this vertical stack.

We complete the proof by noting that

$\begin{array}{ll}{H}_{i,j}& =\sum _k{A}_{i,k}^h\cdot {H^{\mathrm{\prime }\mathrm{\prime }}}_{k,j}=\sum _k\sum _r{A}_{i,k}^r\cdot {H^{\mathrm{\prime }}}_{r,k,j}\\ & =\sum _k\sum _r{A}_{i,k}^r\sum _m{X}_{r,k,m}\cdot W_{r,m,j}=\sum _r\sum _m{W}_{r,m,j}\sum _k{A}_{i,k}^r{X}_{r,k,m}\\ & =\sum _r\sum _m{W}_{r,m,j}{G^{\mathrm{\prime }}}_{r,i,m}=\sum _r{G^{\mathrm{\prime }\mathrm{\prime }}}_{r,i,j}=G_{i,j}.\end{array}$

For the following theorem we define a relational graph $\mathcal{G}$ as a directed graph with labelled nodes and edges. R is the number of relations in the graph and N is the number of nodes in the graph. $\mathcal{N}$ is a set of nodes in graph $\mathcal{G}$. Target nodes, $\mathcal{T}\subseteq \mathcal{N}$, are nodes in the graph with class labels. $h_i^l$ is the vector representation of node i at layer l, and $W_r$ is a relation-specific weight matrix. k is the number of RGCN layers. Non-linear activation function $\sigma$ is an element-wise operation. ${\mathcal{N}}_r(i)$ is the collection of the incoming neighbors of node i with the relation r.

Theorem 2. For a k-layer RGCN model, the input representation $h_m^0$ of a node m more than k hops away from any node in $\mathcal{T}$ is not involved in the computation of $h_i^k$ for any node $i\in \mathcal{T}$.

Proof. We will prove this by induction on k.

Base case (k = 1): Let m be a node which is more than 1 hop away from all nodes in $\mathcal{T}$. We need to show that the computation of $h_i^1$, $i\in \mathcal{T}$ does not involve $h_m^0$.

By definition, $h_i$ is computed as

${\mathbf{h}}_i^1=\sigma \left[\sum _r^R\sum _{j\in {\mathcal{N}}_r(i)}{\displaystyle \frac{1}{|{\mathcal{N}}_r(i)|}}{W}_r{h}_j^0\right].$

The only node representations involved are $h_j^0$ for $j\in {\mathcal{N}}_r(i)$, which are in the 1-hop neighborhood of i.

Inductive step: We assume that the theorem holds for k = l. Let node m be more than l + 1 hops from all nodes in $\mathcal{T}$. We aim to show that for $i\in \mathcal{T}$, $h_m^0$ is not involved in the computation of $h_i$.

The output representation of i is computed as

${\mathbf{h}}_i^{l+1}=\sigma \left[\sum _r^R\sum _{j\in {\mathcal{N}}_r(i)}{\displaystyle \frac{1}{|{\mathcal{N}}_r(i)|}}{W}_r{h}_j^l\right].$

Here j is one hop from i so must be more than l hops from m (if j were l hops from m we could get from $i\to j\to m$ in l + 1 hops). From this we can apply our inductive assumption and conclude that $h_m^0$ is not involved in the computation of $h_j^l$.

We can now see that the only node representations involved in the computation of $h_i^{l+1}$ are $h_j^l$, none of which were computed using $h_m^0$.

Notations

We use lowercase letters to denote scalars (e.g., x), bold lowercase letters for vectors (e.g., w), uppercase letters for matrices (e.g., A), and caligraphic letters for sets (e.g., $\mathcal{G}$). We also use uppercase letters for referring to dimensions and lowercase letters for indexing over those dimensions (e.g., ${\sum }_{b=1}^B$). The dimensions of vectors, matrices and tensors are reported in Table A1.

General experimental setup

All experiments were performed on a single-node machine with an Intel(R) Xeon(R) Gold 5118 (2.30 GHz, 12 cores) CPU and 64 GB of RAM. GPU experiments used a Nvidia GeForce GTX 1080 Ti GPU. We used the Adam optimiser (Kingma & Ba, 2015) with a learning rate of 0.01. The hyperparameters for the reproduction experiements were taken from the original work by Schlichtkrull et al. (2018). However, some hyperparameters were not report in the article, so we obtained the missing configurations from the original codebase and personal communication with the authors. For the new models (c-RGCN & e-RGCN), the hyperparameters were tuned using validation sets. However, we did not perform an extensive systematic hyperparameter search. For reproducibility, we provide an extension description of the hyperparameters that we have used in node classifications and link predictions in YAML files under the configs directory on the project GitHub page: https://github.com/thiviyanT/torch-rgcn.

We ran the original implementation of the link prediction model (https://github.com/MichSchli/RelationPrediction) on the FB15-237 dataset. The exact hyperparameters for WN18 and FB15k experiments were not available in the original codebase. We used Tensorflow 1.4, Theano 1.0.5 and CUDA 8.0.44. This replication required a Nvidia Titan RTX GPU with 24 GB of GPU memory, but model training and inference was performed on the CPU.

Schlichtkrull initialisation

The initialization used in the link prediction models in Schlichtkrull et al. (2018) differs slightly from the more standard Glorot initialization (Glorot & Bengio, 2010).

(17) $$std=\mathrm{\Phi }\times {\displaystyle \frac{3}{\sqrt{\mathrm{f}\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{n}+\mathrm{f}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{u}\mathrm{t}}},}$$

Here, gain ( $\mathrm{\Phi }$) is a constant that is used to scale the standard deviation according the applied non-linearity. std is used to sample random points from either a standard normal or uniform distribution. We refer to this scheme as Schlichtkrull initialisation. When gain is not required, $\mathrm{\Phi }$ is set to 1.0.

FB-Toy link prediction

We also performed link prediction experiment on the FB-Toy dataset. The mean and standard error of the link prediction results are reported in Table A2.