Year: 2021

Abstract:
This paper studies node classification in the inductive setting, i.e., aiming to learn a model on labeled training graphs and generalize it to infer node labels on unlabeled test graphs. This problem has been extensively studied with graph neural networks (GNNs) by learning effective node representations, as well as traditional structured prediction methods for modeling the structured output of node labels, e.g., conditional random fields (CRFs). In this paper, we present a new approach called the Structured Markov Network (SMN), which combines the advantages of both worlds. SMN defines flexible potential functions of CRFs with GNNs. However, learning such a model is nontrivial as it involves optimizing a maximin game with high-cost inference. Inspired by the underlying connection between the joint and marginal distributions defined by Markov networks, we propose to solve an approximate version of the optimization problem as a proxy, which yields a near-optimal solution, making learning more efficient. Extensive experiments on two settings show our approach outperforms many competitive baselines.

text: This paper studies node classification in the inductive setting, i.e., aiming to learn a model on labeled training graphs and generalize it to infer node labels on unlabeled test graphs. This problem has been extensively studied with graph neural networks (GNNs) by learning effective node representations, as well as traditional structured prediction methods for modeling the structured output of node labels, e.g., conditional random fields (CRFs). In this paper, we present a new approach called the Structured Markov Network (SMN), which combines the advantages of both worlds. SMN defines flexible potential functions of CRFs with GNNs. However, learning such a model is nontrivial as it involves optimizing a maximin game with high-cost inference. Inspired by the underlying connection between the joint and marginal distributions defined by Markov networks, we propose to solve an approximate version of the optimization problem as a proxy, which yields a near-optimal solution, making learning more efficient. Extensive experiments on two settings show our approach outperforms many competitive baselines.

heading: 1 INTRODUCTION
text: Graph-structured data (e.g., World Wide Web) are ubiquitous in the real world, covering a variety of applications. This paper studies node classification, a fundamental problem in graph machine learning. Most existing efforts of node classification focus on the transductive setting (Kipf & Welling, 2017; Veličković et al., 2018), i.e., using a small set of labeled nodes in a graph to classify the rest of nodes. In this paper, we study node classification in the inductive setting (Hamilton et al., 2017), which is receiving growing interest in the graph machine learning community. Given a set of training graphs with all the nodes labeled, we aim to learn a classifier to infer the nodes in unlabeled test graphs.
This problem has been recently studied with graph neural networks (GNNs) (Kipf & Welling, 2017; Hamilton et al., 2017; Gilmer et al., 2017; Veličković et al., 2018). GNNs infer the marginal label distribution of each node by learning useful node representations based on node features and edges. Specifically, the representation of each node is iteratively updated by the representations of its neighbors. Once a GNN is learned on training graphs, it can be further applied to test graphs to infer node labels. Owing to the high capacity of nonlinear neural architectures, GNNs achieve impressive results on many datasets. However, one limitation of GNNs is that they ignore the joint dependency of node labels, and thus node labels are predicted separately without modeling structured output.
Indeed, modeling structured output has been widely explored by the literature of structured prediction (BakIr et al., 2007). Structured prediction methods predict node labels collectively, so the label prediction of each node can be improved according to the predicted labels of neighboring nodes. One representative approach is the conditional random field (CRF) (Lafferty et al., 2001). A CRF models the joint distribution of node labels with Markov networks, and thus training CRFs becomes a learning task in graphical models, while predicting node labels corresponds to an inference task. Typically, the potential functions in CRFs are parameterized as log-linear functions, which suffer from low model capacities. One remedy for this is to define potential functions with GNNs (Ma et al., 2018; Qu et al., 2019). However, most of the effective methods for learning CRFs involve a maximin game (Wainwright & Jordan, 2008; Sutton & McCallum, 2012), making learning often hard to converge, especially when GNNs are used to parameterize potential functions. Besides, as learning CRFs requires doing inference on the graphical models, the combined model requires a long run time.
In this paper, we address these challenges by proposing SMN (Structured Markov Network), which is high in capacity, efficient in learning, and able to model the joint dependency of node labels. SMN is inspired by theoretical works in graphical models (Wainwright & Jordan, 2008), which reveal close connections between the joint label distribution and the node/edge marginal label distribution in a Markov network. Based on that, we approximate the original optimization problem with a proxy problem, where the potential functions in CRFs are defined by combining a collection of node/edge pseudomarginal distributions, which are parameterized by GNNs that satisfy a few simple constraints. This proxy problem can be easily solved by maximizing the data likelihood on each node and edge, which yields a near-optimal joint label distribution on training graphs. Once the model is learned, we apply it to test graphs and run loopy belief propagation (Murphy et al., 1999) to infer node labels. Experiments on two settings against both GNNs and CRFs prove the effectiveness of our approach.

heading: 2 RELATED WORK
text: Graph neural networks (GNNs) perform node classification by learning useful node representations. The idea is to iteratively update the representation of each node according to its own features as well as its neighboring nodes. Well-known models include graph convolutional networks (Kipf & Welling, 2017), graph attention networks (Veličković et al., 2018), and message passing neural networks (Gilmer et al., 2017). Most earlier efforts focus on designing GNNs for node classification in the transductive setting (Yang et al., 2016; Kipf & Welling, 2017; Gao & Ji, 2019; Xhonneux et al., 2020), and many recent works move to study the inductive setting (Hamilton et al., 2017; Gao et al., 2018; Chiang et al., 2019; Li et al., 2019; Chen et al., 2020a; Zeng et al., 2020). Because of high model capacity and efficient end-to-end training, GNNs achieve impressive results on inductive node classification. Despite the great success, GNNs only try to model the marginal distribution of each node label and predict node labels separately without considering their joint dependency. In contrast, SMN models joint distributions of node labels by combining GNNs with CRFs, which can thus predicts node labels collectively to improve the result. See Sec. 4.3 for the detailed discussion.
Another type of approach for inductive node classification is structured prediction, which focuses on modeling the dependency of node labels, so that the predicted node labels are more consistent. One representative approach is structured SVM (Tsochantaridis et al., 2005; Finley & Joachims, 2008; Sarawagi & Gupta, 2008), but it lacks a probabilistic interpretation to handle the uncertainty of the prediction. Another representative probabilistic approach is conditional random field (Lafferty et al., 2001; Sutton & McCallum, 2006), which models the distribution of output spaces by using a Markov network. CRFs have been proven effective in many applications, such as POS tagging (Lafferty et al., 2001), shallow parsing (Sha & Pereira, 2003), image labeling (He et al., 2004), and sequence labeling (Lample et al., 2016; Ma & Hovy, 2016; Liu et al., 2018). Nevertheless, the potential functions in CRFs are typically defined as log-linear functions, suffering from low model capacity.
There are also some recent works trying to combine GNNs and CRFs. Some works use GNNs to solve inference problems in graphical models (Dai et al., 2016; Satorras et al., 2019; Zhang et al., 2020; Chen et al., 2020b; Satorras & Welling, 2020). In contrast, our approach uses GNNs to parameterize the potential functions in CRFs, which is in a similar vein to Ma et al. (2018); Qu et al. (2019); Ma et al. (2019; 2021); Wang et al. (2021). Among them, Ma et al. (2018) and Qu et al. (2019) optimize the pseudolikelihood (Besag, 1975) for model learning, and Wang et al. (2021) optimizes a cross-entropy loss on each single node, which can yield poor approximation of the true joint likelihood (Koller & Friedman, 2009; Sutton & McCallum, 2012). Our approach instead solves a proxy problem, which yields a near-optimal solution to the original problem of maximizing likelihood, and thus achieves superior results. For Ma et al. (2019), they model joint dependency of node labels by designing a generative model, which is evaluated in the transductive setting. Different from them, our approach designs a discriminative model based on CRF, which has weaker assumptions on the underlying generative process, and we focus on the inductive setting. For Ma et al. (2021), they study continuous node labels, which is different from our setting of node classification with discrete labels.
Lastly, learning CRFs has been widely explored in graphical models. Many methods perform learning by solving a maximin game as a surrogate for data likelihood maximization (Sutton & McCallum, 2009; 2012; Wiseman & Kim, 2019). Nevertheless, these maximin games are often hard to optimize. To avoid solving the maximin game, we follow Wainwright et al. (2003) and build an approximate optimization problem as a proxy, which is much easier to solve and yields a near-optimal solution.

heading: 3 PRELIMINARY

heading: 3.1 PROBLEM DEFINITION
text: This paper focuses on inductive node classification (Hamilton et al., 2017), where we are given some labeled training graphs and unlabeled test graphs. Each training graph is represented as (y∗V ,xV , E), where xV and y∗V are features and labels of a set of nodes V , and E is a set of edges. For each test graph (xṼ , Ẽ), only node features xṼ and edges Ẽ are observed. Based on that, we aim to solve:
• Learning. On training graphs, learning a probabilistic model to approximate p(yV |xV , E). • Inference. For each test graph, inferring node labels y∗
Ṽ according to the distribution p(yṼ |xṼ , Ẽ).
The problem can be solved by GNNs and CRFs. Next, we introduce the details of these methods.

heading: 3.2 GRAPH NEURAL NETWORKS
text: For inductive node classification, graph neural networks (GNNs) learn node representations to predict marginal label distributions of nodes. GNNs assume all node labels are independent conditioned on node features and edges, so the joint label distribution is factorized into a set of marginals as below:
pθ(yV |xV , E) = ∏ s∈V pθ(ys|xV , E). (1)
Each marginal distribution pθ(ys|xV , E) is modeled as a categorical distribution over label candidates, and the label probabilities are computed by applying a linear softmax classifier to the representation of node s. In general, node representations are learned via the message passing mechanism (Gilmer et al., 2017), which brings high capacity to GNNs. Also, owing to the factorization in Eq. (1), learning and inference can be easily solved in GNNs, where we simply need to compute loss and make prediction on each node separately. However, GNNs approximate only the marginal label distributions of nodes on training graphs, which may generalize badly and result in poor approximation of node marginal label distributions on test graphs. Also, the labels of different nodes are separately predicted according to their own marginal label distributions, yet the joint dependency of node labels is ignored.

heading: 3.3 CONDITIONAL RANDOM FIELDS
text: In the context of the paper, conditional random fields (CRFs) build graphical models for node classification. A popular model is the pair-wise CRF, which formalizes the joint label distribution as:
pθ(yV |xV , E) = 1
Zθ(xV , E) exp{ ∑ s∈V θs(ys,xV , E) + ∑ (s,t)∈E θst(ys, yt,xV , E)} (2)
where Zθ(xV , E) is the partition function. θs(ys,xV , E) and θst(ys, yt,xV , E) are scalar scores contributed by each node s and each edge (s, t). In practice, these θ-functions can be either defined as simple linear functions or complicated GNNs. To make the notation concise, we will omit xV and E in the θ-functions, e.g., simplifying θs(ys,xV , E) as θs(ys). With these θ-functions, CRFs are able to model the joint dependency of node labels and therefore achieve structured prediction.
However, learning CRFs to maximize likelihood pθ(y∗V |xV , E) on training graphs is nontrivial in general, as the partition function Zθ(xV , E) is typically intractable in graphs with loops. Thus, a major line of research instead optimizes a maximin game equivalent to likelihood maximization (Wainwright & Jordan, 2008). The maximin game for each training graph (y∗V ,xV , E) is formalized as follows:
max θ
log pθ(y ∗ V |xV , E) =max θ min q L(θ, q), with L(θ, q) =∑ s∈V {θs(y∗s )− Eqs(ys)[θs(ys)]}+ ∑ (s,t)∈E {θst(y∗s , y∗t )− Eqst(ys,yt)[θst(ys, yt)]} −H[q(yV )]. (3)
Here, q(yV ) is a variational distribution on node labels, qs(ys) and qst(ys, yt) are its marginal distributions on nodes and edges. H[q(yV )] := −Eq(yV )[log q(yV )] is the entropy of q(yV ). Given the maximin game, q and θ can be alternatively optimized via coordinate descent (Sutton & McCallum, 2012). In each iteration, we first update the node and edge marginals {qs(ys)}s∈V , {qst(ys, yt)}(s,t)∈E
towards those defined by pθ. This can be done by MCMC, but the time cost is high, so approximate inference is often used, such as loopy belief propagation (Murphy et al., 1999). After q is optimized, we further update θ-functions with the node and edge marginals defined by q via gradient descent.
The optimal θ-functions are characterized by the following moment-matching conditions:
pθ(ys|xV , E) = Iy∗s {ys} ∀s ∈ V, pθ(ys, yt|xV , E) = I(y∗s ,y∗t ){(ys, yt)} ∀(s, t) ∈ E, (4) where Ia{b} is an indicator function whose value is 1 if a = b and 0 otherwise. See Sec. A and Sec. B in appendix for detailed derivation of the maximin game as well as the moment-matching conditions.
Once the θ-functions are learned, they can be further applied to each test graph (xṼ , Ẽ) to predict the joint label distribution as pθ(yṼ |xṼ , Ẽ). Then the best label assignment y∗Ṽ can be inferred by using approximate inference algorithms, such as loopy belief propagation (Murphy et al., 1999).
The major challenge of CRFs lies in learning. On the one hand, learning relies on inference, meaning that we have to update {qs(ys)}s∈V , {qst(ys, yt)}(s,t)∈E to approximate the node and edge marginals of pθ at each step, which can be expensive. On the other hand, as learning involves a maximin game and the optimal q of the inner minimization problem in Eq. (3) is intractable, we can only maximize an upper bound of the likelihood function for θ, making learning unstable. The problem becomes even more severe when θ is parameterized by highly nonlinear neural models, e.g. GNNs.

heading: 4 MODEL
text: Next, we introduce SMN. SMN combines GNNs and CRFs by parameterizing potential functions in CRFs with GNNs, and thus enjoys high capacity and can model the joint dependency of node labels.
However, as elaborated in Sec. 3.3, learning such a model on training graphs is challenging due to the maximin game in optimization. Inspired by the connection between the joint and marginal distributions of CRFs, we instead construct a new optimization problem, which serves as a proxy for model learning. Compared with the original maximin game, the proxy problem is much easier to solve, where we can simply train two GNNs to approximate the marginal label distributions on nodes and edges, and further combine these pseudomarginals (defined in Prop. 1) into a near-optimal joint label distribution. This joint label distribution can be further refined by optimizing the maximin game, although it is optional and often unnecessary, as this distribution is often close enough to the optimal one. With this proxy problem for model learning, learning becomes more stable and efficient.
Afterwards, the learned model is used to predict the joint label distribution on test graphs. Then we run loopy belief propagation to infer node labels. Now, we introduce the details of our approach.

heading: 4.1 LEARNING
text: The learning task aims at training θ to maximize the log-likelihood function log pθ(y∗V |xV , E) for each training graph (y∗V ,xV , E), which is highly challenging. Therefore, instead of directly optimizing this goal, we solve an approximate version of the problem as a proxy, which is training a node GNN and an edge GNN to maximize the log-likelihood of observed labels on nodes and edges.
The Proxy Problem. The proxy problem is inspired by Wainwright & Jordan (2008), which points out that the marginal label distributions on nodes and edges defined by a Markov network have inherent connections with the joint distribution. This connection is stated in the proposition below.
Proposition 1 Consider a set of nonzero pseudomarginals {τs(ys)}s∈V and {τst(ys, yt)}(st)∈E which satisfy ∑ ys τst(ys, yt) = τt(yt) and ∑ yt τst(ys, yt) = τs(ys) for all (s, t) ∈ E.
If we parameterize the θ-functions of pθ in Eq. (2) in the following way:
θs(ys) = log τs(ys) ∀s ∈ V, θst(ys, yt) = log τst(ys, yt)
τs(ys)τt(yt) ∀(s, t) ∈ E, (5)
then {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E are specified by a fixed point of the sum-product loopy belief propagation algorithm when applied to the joint distribution pθ, which implies that:
τs(ys) ≈ pθ(ys) ∀s ∈ V, τst(ys, yt) ≈ pθ(ys, yt) ∀(s, t) ∈ E. (6)
The proof is provided in Sec. C. With the proposition, we observe that if we parameterize the θfunctions by combining a set of pseudomarginals {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E in the way defined by Eq. (5), then those pseudomarginals can well approximate the true marginals of the joint distribution pθ, i.e., τs(ys) ≈ pθ(ys) and τst(ys, yt) ≈ pθ(ys, yt) for all nodes s and edges (s, t). Given this precondition, if we further have τs(ys) ≈ Iy∗s {ys} and τst(ys, yt) ≈ I(y∗s ,y∗t ){(ys, yt)}, then the moment-matching conditions in Eq. (4) for the optimal θ-functions are roughly satisfied. This implies the joint distribution pθ(yV |xV , E) derived in this way is a near-optimal one. With the observation, rather than directly using GNNs to parameterize the θ-functions, we use a node GNN and an edge GNN to parameterize the pseudomarginals {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E . For the pseudomarginal τs(ys) on node s, we apply the node GNN to node features xV and edges E, yielding a representation us for node s. Then we apply a softmax classifier to us to compute τs(ys):
{us}u∈V = GNNnode(xV , E), τs(ys) = softmax(f(us))[ys], (7)
where f maps a node representation to a |Y|-dimensional logit and Y is the node label set. Similarly, we apply the edge GNN to compute a representation vs for each node s, and model τst(ys, yt) as:
{vs}s∈V = GNNedge(xV , E) τst(ys, yt) = softmax(g(vs,vt))[ys, yt], (8)
where g is a function mapping a pair of representations to a (|Y| × |Y|)-dimensional logit. Given the parameterization, we construct the following problem as a proxy for learning θ-functions:
min τ,θ ∑ s∈V d ( Iy∗s {ys}, τs(ys) ) + ∑ (s,t)∈E d ( I(y∗s ,y∗t ){(ys, yt)}, τst(ys, yt) ) ,
subject to θs = log τs(ys), θst(ys, yt) = log τst(ys, yt)
τs(ys)τt(yt) , and ∑ ys τst(ys, yt) = τt(yt), ∑ yt τst(ys, yt) = τs(ys),
(9)
for all nodes and edges, where d can be any divergence measure between two distributions. By solving the above problem, {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E will be valid pseudomarginals which can well approximate the true labels, i.e., τs(ys) ≈ Iy∗s {ys} and τst(ys, yt) ≈ I(y∗s ,y∗t ){(ys, yt)}. Then according to the constraint in the second line of Eq. (9), θ-functions are formed in a way to enable τs(ys) ≈ pθ(ys) and τst(ys, yt) ≈ pθ(ys, yt) as stated in the Prop. 1. Combining these two sets of formula results in pθ(ys) ≈ Iy∗s {ys} and pθ(ys, yt) ≈ Iy∗s {ys}. We see that the moment-matching conditions in Eq. (4) for the optimal joint label distribution are roughly achieved, implying that the derived joint distribution pθ(yV |xV , E) is a near-optimal solution to the original learning problem. One good property of the proxy problem is that it can be solved easily. The last consistency constraint (i.e. ∑ ys τst(ys, yt) = τt(yt) and ∑ yt τst(ys, yt) = τs(ys)) can be ignored during optimization, since by optimizing the objective function, the optimal pseudomarginals τ should well approximate the observed node and edge marginals, i.e., τs(ys) ≈ Iy∗s {ys} and τst(ys, yt) ≈ I(y∗s ,y∗t ){(ys, yt)}, and hence τ will almost naturally satisfy the consistency constraint. We also tried some constrained
optimization methods to handle the consistency constraint, but they yield no improvement. See Sec. D of appendix for more details. Thus, we can simply train the pseudomarginals parameterized by GNNs to approximate the true node and edge labels on training graphs, i.e., minimizing d(Iy∗s {ys}, τs(ys)) and d(I(y∗s ,y∗t ){(ys, yt)}, τst(ys, yt)). Then we build θ-functions as in Eq. (5) to obtain a near-optimal joint distribution. In practice, we choose d to be the KL divergence, yielding an objective for τ as:
max τ ∑ s∈V log τs(y ∗ s ) + ∑ (s,t)∈E log τst(y ∗ s , y ∗ t ). (10)
This objective function is very intuitive, where we simply try to optimize the node GNN and edge GNN to maximize the log-likelihood function of the observed labels on nodes and edges.
Refinement. By solving the proxy problem, we can obtain a near-optimal joint distribution. In practice, we observe that when we have a large amount of training data, further refining this joint distribution by solving the maximin game in Eq. (3) for a few iterations can lead to further improvement. Formally, each iteration of refinement has two steps. In the first step, we run sum-product loopy belief propagation (Murphy et al., 1999), which yields a collection of node and edge marginals (i.e., {qs(ys)}s∈V and {qst(ys, yt)}(s,t)∈E) as approximation to the marginals defined by pθ. In the second step, we update the θ-functions parameterized by the node and edge GNNs to maximize:∑
s∈V
{ θs(y ∗ s )− Eqs(ys)[θs(ys)] } + ∑ (s,t)∈E { θst(y ∗ s , y ∗ t )− Eqst(ys,yt)[θst(ys, yt)] } . (11)
Intuitively, we treat the true label y∗s and (y ∗ s , y ∗ t ) of each node and edge as positive examples, and encourage the θ-functions to raise up their scores. Meanwhile, those labels sampled from qs(ys) and qst(ys, yt) act as negative examples, and the θ-functions are updated to decrease their scores.

heading: 4.2 INFERENCE
text: After learning, we apply the node and edge GNNs to each test graph (xṼ , Ẽ) to compute the θfunctions, which are integrated into an approximate joint label distribution pθ(yṼ |xṼ , Ẽ). Then we use this distribution to infer the best label y∗s̃ for each node s̃ ∈ Ṽ , where two settings are considered. Node-level Accuracy. Typically, we care about the node-level accuracy, i.e., how likely we can correctly classify a node in test graphs. Intuitively, the best label y∗s̃ for each test node s̃ ∈ Ṽ should be predicted as y∗s̃ = argmaxys̃ pθ(ys̃|xṼ , Ẽ), where pθ(ys̃|xṼ , Ẽ) is the marginal label distribution of node s̃ induced by the joint pθ(yṼ |xṼ , Ẽ). In practice, the exact marginal is intractable, so we apply loopy belief propagation (Murphy et al., 1999) for approximate inference. For each edge (s̃, t̃) in test graphs, we introduce a message function mt̃→s̃(ys̃) and iteratively update all messages as:
mt̃→s̃(ys̃) ∝ ∑ yt̃ {exp(θt̃(yt̃) + θs̃t̃(ys̃, yt̃)) ∏ s̃′∈N(t̃)\s̃ ms̃′→t̃(yt̃)}, (12)
where N(s̃) denotes the set of neighboring nodes for node s̃. Once the above process converges or after sufficient iterations, the label of each node s̃ can be inferred in the following way:
y∗s̃ = argmax ys̃
[exp(θs̃(ys̃)) ∏
t̃∈N(s̃)
mt̃→s̃(ys̃)]. (13)
Graph-level Accuracy. In some other cases, we might care about the graph-level accuracy, i.e., how likely we can correctly classify all nodes in a given test graph. In this case, the best prediction of node labels is given by y∗
Ṽ = argmaxyṼ p(yṼ |xṼ , Ẽ). This problem can be approximately solved by the
max-product variant of loopy belief propagation, which simply replaces the sum over yt̃ in Eq. (12) with max (Weiss & Freeman, 2001). Afterwards, the best node label can be still decoded via Eq. (13).

heading: 4.3 DISCUSSION
text: Next, we discuss the advantages of SMN over existing GNN methods. (1) Compared with GNNs which only learn marginal label distributions on nodes, SMN approximates marginals on both nodes and edges. These node and edge marginals mutually benefit each other, allowing SMN to generalize better to test graphs. (2) Besides, SMN tries to model the joint distribution of node labels, which allows for collective inference to improve node label prediction. It also allows SMN to choose proper inference algorithms for different purposes, e.g., maximizing node-level or graph-level accuracy.

heading: 5 EXPERIMENT

heading: 5.1 DATASETS
text: We consider datasets in two settings, which focus on node-level and graph-level accuracy respectively.
Node-level Accuracy. The node-level accuracy measures how likely a model can predict the correct label of a node in test graphs. We use the PPI dataset (Zitnik & Leskovec, 2017; Hamilton et al., 2017), which has 20 training graphs. To make the dataset more challenging, we also try using only the first 1/2/10 training graphs, yielding another three datasets PPI-1, PPI-2, and PPI-10. Besides, we also build a DBLP dataset from the citation network 1 in Tang et al. (2008). Papers from eight conferences are treated as nodes, and we split them into three categories for classification according to conference domains 2. For each paper, we compute the mean GloVe embedding 3 (Pennington et al., 2014) of words in the title and abstract as node features. The training/validation/test graph is formed as the citation graph of papers published before 1999, from 2000 to 2009, after 2010 respectively.
Graph-level Accuracy. The graph-level accuracy measures how likely a model can correctly classify all the nodes for a given test graph. We construct three datasets from the Cora, Citeseer, and Pubmed datasets used for transductive node classification (Yang et al., 2016). Each raw dataset has a single graph. For each training/validation/test node of the raw dataset, we treat its ego network 4 as a training/validation/test graph. We denote the datasets as Cora*, Citeseer*, Pubmed*.

heading: 5.2 COMPARED ALGORITHMS
text: Graph Neural Networks. For GNNs, we choose a few well-known model architectures for comparison, including GCN (Kipf & Welling, 2017), GraphSage (Hamilton et al., 2017), GAT (Veličković et al., 2018), Graph U-Net (Gao & Ji, 2019) and GCNII (Chen et al., 2020a).
Conditional Random Fields. For CRFs, we consider three variants. (1) CRF-linear. This variant uses linear θ-functions in Eq. (2), which takes the features on nodes and edges for computation. (2) CRF-GNN. This variant parameterizes the θ-functions as θs(ys) = f(us) and θst(ys, yt) = g(vs,vt), with f and g defined in Eq. (7) and Eq. (8), where the node representations are generated by different GNN architectures (e.g., CRF-GAT). We train these models via the maximin game as in Eq. (3) with sum-product loopy belief propagation. (3) GMNN. We also consider GMNN (Qu et al., 2019), an approach combining GNNs and CRFs, which optimizes the pseudolikelihood function for learning.
Our Approach. For SMNs, we try different GNN architectures for defining the node and edge GNNs (e.g., SMN-GAT). By default, we only solve the proxy problem without performing refinement. We systematically compare the results with and without refinement in part 2 of Sec. 5.5.

heading: 5.3 EVALUATION METRICS
text: On Cora*, Citeseer*, and Pubmed*, we report the percentage of test graphs where all the nodes are correctly classified (i.e., graph-level accuracy). On DBLP and PPI, we report accuracy and micro-F1 based on the percentage of test nodes which are correctly classified (i.e., node-level accuracy). For Cora*, Citeseer*, and Pubmed*, we run each compared method with 10 different seeds to report the mean accuracy and the standard deviation. For DBLP and PPI, we run each method with 5 seeds.

heading: 5.4 EXPERIMENTAL SETUP
text: For GNNs, by default we use the same architectures (e.g., number of neurons, number of layers) as used in the original papers. Adam (Kingma & Ba, 2015) is used for training. For the edge GNN in Eq. (8), we add a hyperparameter γ to control the annealing temperature of the logit g(vs,vt) before the softmax function during belief propagation. Empirically, we find that max-product belief propagation works better than the sum-product variant in most cases, so we use the max-product version by default. By default, we do not run refinement when training SMNs. See Sec. F for details.
1https://originalstatic.aminer.cn/misc/dblp.v12.7z 2ML: ICML/NeurIPS. CV: ICCV/CVPR/ECCV. NLP: ACL/EMNLP/NAACL. 3http://nlp.stanford.edu/data/glove.6B.zip 4The local subgraph formed by a node and its direct neighbors.

heading: 5.5 RESULTS
text: 1. Comparison with other methods. The main results in the two settings are presented in Tab. 1 and Tab. 2. Compared against different GNN models, our approach achieves consistent improvement (the relative underperformance of SMN-GCN and SMN-SAGE is related to the capacity of the backbone GNNs and is explained in Sec. G.1) by using these GNNs as backbone networks for approximating marginal label distributions on nodes and edges, which demonstrates SMNs are able to model the structured output of node labels by combining with CRFs, and thus achieve better results.
Besides, SMNs also achieve superior results to CRF-GNNs which are trained by directly solving the maximin game in Eq. (3), as well as GMNN which optimizes the pseudolikelihood function. This observation proves the advantage of our proposed proxy optimization problem for learning CRFs.
Table 3: Run time comparison (in sec).
Algorithm DBLP PPI GAT 23.15 460.81 CRF (GAT) 500.43 27136.90 SMN(GAT) 46.86 962.92
Table 4: Micro-F1 with and w/o refinement (in %).
Algorithm Refine PPI-2 PPI-10 PPI SMN- w/o 71.52 ± 0.21 94.41 ± 0.21 98.38 ± 0.10 GAT with 71.58 ± 0.20 94.63 ± 0.20 98.68 ± 0.09
SMN- w/o 73.93 ± 0.08 91.99 ± 0.04 97.56 ± 0.03 GraphSAGE with 73.68 ± 0.10 92.49 ± 0.02 97.77 ± 0.02
2. Effect of refinement. By solving the proxy optimization problem in Eq. (9), we can obtain a near-optimal joint label distribution on training graphs, based on which we may optionally refine the distribution with the maximin game in Eq. (3). Next, we study the effect of refinement, and we present the results in Tab. 4. By only solving the proxy problem, our approach already achieves impressive results, showing that the proxy problem can well approximate the original learning problem. Only on datasets with sufficient labeled data (e.g., PPI-10, PPI), refinement leads to some improvement.
3. Model architecture. SMN uses a node GNN and an edge GNN for computing node and edge marginals independently. In practice, we can also use a shared GNN for both node and edge marginals. We show results of this variant in Tab. 6, where it also achieves significant improvement over GNNs.
4. Efficiency comparison. We have seen SMNs achieve better classification results than GNNs and CRFs. Next, we further compare their efficiency by showing the run time on DBLP and PPI.
Table 5: Comparison of learning methods (in %).
Algorithm Cora* Citeseer* PPI-10
Maximin Game 49.10 ± 3.80 42.89 ± 1.30 54.98 ± 1.13
Pseudolikelihood 54.30 ± 1.15 48.46 ± 1.06 90.72 ± 5.28
Proxy Problem 58.78 ± 1.21 49.02 ± 0.78 95.87 ± 0.02
Table 6: Micro-F1 of model variants (in %).
Algorithm PPI-1 PPI-2 PPI-10
GAT 60.72 ± 0.25 68.55 ± 0.30 93.53 ± 0.24 SMN-GAT 64.02 ± 0.40 72.37 ± 0.18 94.41 ± 0.21node and edge GNNs SMN-GAT 63.72 ± 0.38 70.99 ± 0.25 95.19 ± 0.15a shared GNN
For PPI, which has 121 labels, we only report the training times on a single label. We use GAT as the backbone network for CRFs and SMNs. GAT and CRF are trained for 1000 epochs to ensure convergence. For the SMN, we train the node GNN and edge GNN for node/edge classification as in Eq. (10) with 1000 epochs. The run times are presented in Tab. 3. SMNs take twice as long for training than GAT, as a SMN needs to train a node GNN and an edge GNN. Compared with CRFs, SMNs are much more efficient. The reason is that the proxy optimization problem is much easier to solve, as it does not requires running expensive approximate inference algorithms.
5. Comparison of learning methods. Next, we investigate different methods for learning SMNs, including directly solving the maximin game, optimizing pseudolikelihood, and solving our proposed proxy problem. We show the results for optimizing SMN-GAT in Tab. 5. We see solving maximin game yields poor results due to unstable training. Although the pseudolikelihood method performs much better, the result is still unsatisfactory as it is not a good approximation of the true likelihood. By solving our proposed proxy problem, SMN achieves the best result, which proves its effectiveness.
6. Convergence analysis. To better illustrate the advantage of the proxy problem for learning CRFs, we look into the training curves of SMNs, SMNs w/o proxy, and CRFs when optimizing the maximin game in Eq. (3). For SMNs, we optimize the node and edge GNNs on the proxy optimization problem in Eq. (9) before doing refinement with the maximin game, while for SMNs w/o proxy we directly perform refinement with the maximin game without solving the proxy problem. We show the results in Fig. 2. CRFs and SMNs w/o proxy suffer from high variance and low accuracy. In contrast, owing to the near-optimal joint distribution found by solving the proxy problem, SMNs get much higher accuracy with lower variance even without refinement (see initial results of SMNs at epoch 0). Also, the refinement process quickly converges after only a few epochs, showing good efficiency of SMNs.
7. Case study. To intuitively see how SMNs outperform GNNs, we conduct some case studies on Cora*. We use GAT as backbone networks, and show the prediction made by the GAT (the node GNN), the edge GNN, and SMN in Fig. 3. In all three cases shown in the figure, GAT (left column) makes inconsistent predictions on linked nodes, as it fails to model the structured output. The edge GNN (middle column) also makes a mistake in the bottom case. Finally, by combining GAT and edge GNN with a CRF, the SMN (right column) is able to predict the correct labels for all nodes.

heading: 6 CONCLUSION
text: This paper studied inductive node classification, and we proposed SMN to combine GNNs and CRFs. Inspired by the connection of joint and marginal distributions defined by Markov networks, we designed a proxy optimization problem for efficient model learning. In the future, we plan to explore techniques on training GANs (Goodfellow et al., 2014) to better optimize the maximin game for learning. In addition, SMNs model joint dependency of node labels by defining potential functions on nodes and edges, and we also plan to further explore high-order local structures, e.g., triangles.

heading: A DERIVATION OF THE MAXIMIN GAME
text: As discussed in the Sec. 3, optimizing the joint label distribution pθ to maximize the log-likelihood log pθ(y ∗ V |xV , E) on a training graph (y∗V ,xV , E) is equivalent to solving a maximin game. In this section, we provide the detailed derivation.
Let ψθ(yV ,xV , E) be the potential function as below:
ψθ(yV ,xV , E) = exp ∑ s∈V θs(ys,xV , E) + ∑ (s,t)∈E θst(ys, yt,xV , E)  . (14) For each training graph (y∗V ,xV , E), we aim at maximizing the following log-likelihood function:
log pθ(y ∗ V |xV , E) = log
1
Zθ(xV , E) ψθ(y
∗ V ,xV , E)
= logψθ(y ∗ V ,xV , E)− logZθ(xV , E)
= logψθ(y ∗ V ,xV , E)− log ∑ yV ψθ(yV ,xV , E).
(15)
However, the term log ∑
yV ψθ(yV ,xV , E) is computationally intractable, as we need to sum over
all the possible yV . To solve the problem, we introduce a variational joint distribution q(yV ) defined on all node labels yV , and use the Jensen’s inequality to derive an estimation of the term log ∑
yV ψθ(yV ,xV , E) as follows:
log ∑ yV ψθ(yV ,xV , E) = logEq(yV ) [ ψθ(yV ,xV , E) q(yV ) ] ≥ Eq(yV ) [ log ψθ(yV ,xV , E)
q(yV ) ] = Eq(yV )[logψθ(yV ,xV , E)]− Eq(yV )[log q(yV )].
(16)
It is easy to see that the equation holds if and only if q(yV ) = pθ(yV |xV , E), and hence:
log ∑ yV ψθ(yV ,xV , E) = max q(yV ) { Eq(yV )[logψθ(yV ,xV , E)]− Eq(yV )[log q(yV )] } . (17)
By taking the above result into Eq. (15), we obtain:
log pθ(y ∗ V |xV , E) = logψθ(y∗V ,xV , E)− log ∑ yV ψθ(yV ,xV , E)
= min q(yV )
{ logψθ(y ∗ V ,xV , E)− Eq(yV )[logψθ(yV ,xV , E)] + Eq(yV )[log q(yV )] } .
(18)
As ψθ(yV ,xV , E) = exp{ ∑ s∈V θs(ys,xV , E) + ∑ (s,t)∈E θst(ys, yt,xV , E)}, we have:
log pθ(y ∗ V |xV , E) = min q L(θ, q), (19)
with:
L(θ, q) = −H[q(yV )] + ∑
(s,t)∈E {θst(y∗s , y∗t )− Eqst(ys,yt)[θst(ys, yt)]}+ ∑ s∈V {θs(y∗s )− Eqs(ys)[θs(ys)]}. (20)
Therefore, optimizing θ to maximize the log-likelihood function is equivalent to solving the following maximin game:
max θ
log pθ(y ∗ V |xV , E) = max θ min q L(θ, q). (21)

heading: B DERIVATION OF THE MOMENT-MATCHING CONDITIONS
text: In the CRF model defined in the preliminary section, the parameter θ consists of the output values of all θ-functions. In other words, θ = {θs(ys)}ys∈Y,s∈V ∪ {θst(ys, yt)}ys∈Y,yt∈Y,s∈V , where Y is the set of all the possible node labels.
By definition, pθ(yV |xV , E) belongs to the exponential family. According to the property of exponential family distributions, log pθ(y∗V |xV , E) is strictly concave with respect to θ. Therefore, the optimal θ is unique, which is characterized by the condition of ∂∂θ log pθ(y ∗ V |xV , E) = 0. Formally, ∂∂θ log pθ(y ∗ V |xV , E) can be computed as below:
∂
∂θs(ŷs) log pθ(y
∗ V |xV , E) =
∂
∂θ logψθ(y
∗ V ,xV , E)−
∂
∂θ logZθ(xV , E). (22)
For ∂∂θ logZθ(xV , E), we have: ∂
∂θ logZθ(xV , E) =
∂ ∂θ log ∑ yV ψθ(yV ,xV , E)
=
∑ yV ∂ ∂θψθ(yV ,xV , E)∑
yV ψθ(yV ,xV , E)
=
∑ yV ψθ(yV ,xV , E) ∂ ∂θ logψθ(yV ,xV , E)∑
yV ψθ(yV ,xV , E)
= ∑ yV [ ψθ(yV ,xV , E)∑ y′V ψθ(y′V ,xV , E) ∂ ∂θ logψθ(yV ,xV , E) ]
= ∑ yV [ ψθ(yV ,xV , E) Zθ ∂ ∂θ logψθ(yV ,xV , E) ] = Epθ(yV |xV ,E) [ ∂
∂θ logψθ(yV ,xV , E)
] .
(23)
By combining the above two equations, we have: ∂
∂θ log pθ(y
∗ V |xV , E) =
∂
∂θ logψθ(y
∗ V ,xV , E)− Epθ(yV |xV ,E)
[ ∂
∂θ logψθ(yV ,xV , E)
] . (24)
The potential function ψθ above is defined as ψθ(yV ,xV , E) = exp{ ∑ s∈V θs(ys,xV , E) +∑
(s,t)∈E θst(ys, yt,xV , E)}. If we consider each specific scalar θs(ŷs), and taking the derivative with respect to the scalar to 0, we obtain:
0 = ∂
∂θs(ŷs) log pθ(y
∗ V |xV , E)
= ∂
∂θs(ŷs) logψθ(y
∗ V ,xV , E)− Epθ(yV |xV ,E)
[ ∂
∂θs(ŷs) logψθ(yV ,xV , E) ] = Iy∗s {ŷs} [ ∂
∂θs(ŷs) θs(ŷs) ] − ∑ yV pθ(yV |xV , E) [ Iy∗s {ŷs} ∂ ∂θs(ŷs) θs(ŷs) ] = Iy∗s {ŷs} [ ∂
∂θs(ŷs) θs(ŷs)
] − pθ(ŷs|xV , E) [ ∂
∂θs(ŷs) θs(ŷs) ] = Iy∗s {ŷs} − pθ(ŷs|xV , E),
(25)
which implies pθ(ŷs|xV , E) = Iy∗s {ŷs} for the optimal θ. Moreover, this equation holds for all s ∈ V and all ŷs ∈ Y . Similarly, for each scalar θst(ŷs, ŷt), we have that ∂∂θst(ŷs,ŷt) log pθ(y ∗ V |xV , E) = 0 is equivalent to pθ(ŷs, ŷt|xV , E) = Iy∗s ,y∗t {ŷs, ŷt}. This equation holds for all (s, t) ∈ E and all the choices of (ŷs, ŷt) ∈ Y × Y . Therefore, the optimal θ-functions are characterized by the moment-matching conditions as below:
pθ(ys|xV , E) = Iy∗s {ys} ∀s ∈ V, pθ(ys, yt|xV , E) = Iy∗s ,y∗t {ys, yt} ∀(s, t) ∈ E. (26)

heading: C PROOF OF PROPOSITION 1
text: Next, we prove Prop. 1. We first restate the proposition as follows:
Proposition For a node set V and an edge set E, consider a collection of nonzero pseudomarginals {τs(ys)}s∈V and {τst(ys, yt)}(st)∈E which satisfy ∑ ys τst(ys, yt) = τt(yt) and ∑ yt τst(ys, yt) = τs(ys) for all (s, t) ∈ E. If we parameterize the θ-functions in pθ as follows:
θs(ys) = log τs(ys) ∀s ∈ V, θst(ys, yt) = log τst(ys, yt)
τs(ys)τt(yt) ∀(s, t) ∈ E, (27)
then {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E are specified by a fixed point of the sum-product loopy belief propagation algorithm when applied to the joint distribution pθ, implying:
τs(ys) ≈ pθ(ys) ∀s ∈ V, τst(ys, yt) ≈ pθ(ys, yt) ∀(s, t) ∈ E. (28) Proof: To prove the proposition, we first summarize the workflow of the sum-product loopy belief propagation algorithm. In sum-product loopy belief propagation, we introduce a message function mt→s(ys) for each edge (s, t) ∈ E. Then we iteratively update all message functions as follows:
mt→s(ys) ∝ ∑ yt exp(θt(yt) + θst(ys, yt)) ∏ s′∈N(t)\s ms′→t(yt)  , (29) where N(t) represents the set of neighbors for node t.
Once the process converges or after sufficient iterations, the approximation of the node marginals and the edge marginals (i.e., {qs(ys)}s∈V and {qst(ys, yt)}(s,t)∈E) can be recovered by the message functions {mt→s(ys)}(s,t)∈E as follows:
qs(ys) ∝ exp(θs(ys)) ∏
t∈N(s)
mt→s(ys), (30)
qst(ys, yt) ∝ exp(θs(ys) + θt(yt) + θst(ys, yt)) ∏
t′∈N(s)\t
mt′→s(ys) ∏
s′∈N(t)\s
ms′→t(yt). (31)
Next, let us move back to our case, where we parameterize the θ-functions with a set of pseudomarginals as in Eq. (27). For such a specific parameterization of the θ-functions, we claim that one fixed point of Eq. (29) is achieved when mt→s(ys) = 1 for all (s, t) ∈ E. To prove that, we notice that when all the message functions equal to 1, the left side of Eq. (29) is apparently 1. The right side of Eq. (29) can be computed as below:∑
yt
exp(θt(yt) + θst(ys, yt)) ∏
s′∈N(t)\s ms′→t(yt) = ∑ yt exp(θt(yt) + θst(ys, yt))
= ∑ yt exp ( log τt(yt) + log τst(ys, yt) τs(ys)τt(yt) ) = ∑ yt exp ( log τst(ys, yt) τs(ys)
) = ∑ yt τst(ys, yt) τs(ys)
= τs(ys)
τs(ys)
= 1. (32)
We can see that both the left side and the right side of Eq. (29) are 1, and hence {mt→s(ys) = 1}(s,t)∈E specifies a fixed point of sum-product loopy belief propagation. For this fixed point, qs(ys) can be computed as follows:
qs(ys) ∝ exp(θs(ys)) ∏
t∈N(s)
mt→s(ys) = exp(θs(ys)) = τs(ys). (33)
Similarly, we can compute qst(ys, yt) as: qst(ys, yt) ∝ exp(θs(ys) + θt(yt) + θst(ys, yt)) ∏
t′∈N(s)\t
mt′→s(ys) ∏
s′∈N(t)\s
ms′→t(yt)
= exp(θs(ys) + θt(yt) + θst(ys, yt))
= exp ( log τs(ys) + log τt(yt) + log τst(ys, yt)
τs(ys)τt(yt) ) = τst(ys, yt).
(34)
From the above two equations, we can see that {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E are specified by a fixed point (i.e., mt→s(ys) = 1 for all (s, t) ∈ E) of sum-product loopy belief propagation. As sum-product loopy belief propagation often works well in practice to approximate the marginal distributions on nodes and edges, we thus have τs(ys) ≈ pθ(ys) for each node and τst(ys, yt) ≈ pθ(ys, yt) for each edge.

heading: D SOLVING THE PROXY PROBLEM WITH CONSTRAINED OPTIMIZATION
text: The key innovation of our proposed approach is on the proxy optimization problem which is used to approximate the original learning problem. Formally, the proxy optimization problem is stated as:
min τ,θ ∑ s∈V d ( Iy∗s {ys}, τs(ys) ) + ∑ (s,t)∈E d ( I(y∗s ,y∗t ){(ys, yt)}, τst(ys, yt) ) ,
subject to θs = log τs(ys), θst(ys, yt) = log τst(ys, yt)
τs(ys)τt(yt) , and ∑ ys τst(ys, yt) = τt(yt), ∑ yt τst(ys, yt) = τs(ys),
(35)
for all nodes and edges, where d can be any divergence measure between two distributions. In our implementation, we ignore these consistency constraints, i.e., ∑ ys τst(ys, yt) = τt(yt) and∑
yt τst(ys, yt) = τs(ys). This ie because by by optimizing the objective, the obtained pseudomarginals τ would well approximate the observed node and edge marginals, i.e., τs(ys) ≈ Iy∗s {ys} and τst(ys, yt) ≈ I(y∗s ,y∗t ){(ys, yt)}, and hence τ would almost naturally satisfy the constraints. To demonstrate ignoring the consistency constraint makes sense, we also tried a constrained optimization method for solving the proxy problem. Specifically, we add a quadratic term to penalize the inconsistency between τst(ys, yt) and τs(ys) as well as τt(yt), resulting in the following problem:
min τ,θ ∑ s∈V d ( Iy∗s {ys}, τs(ys) ) + ∑ (s,t)∈E d ( I(y∗s ,y∗t ){(ys, yt)}, τst(ys, yt) ) + α
∑ (s,t)∈E [∑ yt { ∑ ys τst(ys, yt)− τt(yt)}2 + ∑ ys { ∑ yt τst(ys, yt)− τs(ys)}2 ] ,
subject to θs = log τs(ys), θst(ys, yt) = log τst(ys, yt)
τs(ys)τt(yt) ,
(36)
for all nodes and edges. Again, d is a divergence measure between two distributions, and we choose to use the KL divergence. α is a hyperparameter deciding the weight of the penalty term.
We conduct empirical comparison of this constrained optimization method and our default implementation where the consistency constraint is ignored. The results are presented in Tab. 7. We can see that the constrained optimization method does not lead to improvement, which shows that ignoring the consistency constraint is empirically reasonable.

heading: E UNDERSTANDING SMNS AS OPTIMIZING A SURROGATE FOR THE LOG-LIKELIHOOD FUNCTION
text: In the model section, we motivate SMNs from the moment-matching conditions of the optimal θfunctions. Specifically, we initialize the θ-functions at a state where the moment-matching conditions are approximately satisfied, yielding a near-optimal joint distribution. Then we further tune the θ-functions to solve the maximin game. Besides this perspective, SMNs can also be understood as optimizing a surrogate for the log-likelihood function. Next, we introduce the details.
Remember that maximizing the log-likelihood function is equivalent to solving a maximin game as:
max θ
log pθ(y ∗ V |xV , E) = max θ min q L(θ, q), with L(θ, q) = −H[q(yV )]
+ ∑ s∈V {θs(y∗s )− Eqs(ys)[θs(ys)]}+ ∑ (s,t)∈E {θst(y∗s , y∗t )− Eqst(ys,yt)[θst(ys, yt)]}. (37)
Here, q(yV ) is a joint distribution on all the node labels. qs(ys) and qst(ys, yt) are the corresponding marginal distributions.
Although the above maximin game is equivalent to the original problem of maximizing likelihood, solving the maximin game is nontrivial. In particular, there are two key challenges, i.e., (1) how to specify constraints to characterize a valid joint distribution q(yV ) and (2) how to compute its entropy H(q) = −Eq(yV )[log q(yV )]. To deal with the challenge, a common practice used in loopy belief propagation is to make the following two approximations:
(1) Instead of specifying constraints to let q(yV ) be a valid joint distribution, we introduce a set of pseudomarginals as approximation to a valid joint distribution. Specifically, these pseudomarginals are denoted as q̃ = {qs(ys)}s∈V ∪ {qst(ys, yt)}(s,t)∈E , and they satisfy ∑ ys qst(ys, yt) = qt(yt)
and ∑ yt qst(ys, yt) = qs(ys) for all (s, t) ∈ E.
(2) We approximate the entropy H(q) with Bethe entropy approximation HBethe(q̃), which is defined as follows:
HBethe(q̃) = − ∑ s∈V Eqs(ys)[log qs(ys)]− ∑ (s,t)∈E Eqst(ys,yt) [ log qst(ys, yt) qs(ys)qt(yt) ] . (38)
With the two approximations, we get the following maximin game as a surrogate for the likelihood maximization problem:
max θ
log pθ(y ∗ V |xV , E) ≈ max θ min q̃ LBethe(θ, q̃), (39)
with: LBethe(θ,q̃) = −HBethe(q̃)
+ ∑
(s,t)∈E {θs,t(y∗s , y∗t )− Eqst(ys,yt)[θs,t(ys, yt)]}+ ∑ s∈V {θs(y∗s )− Eqs(ys)[θs(ys)]}. (40)
This problem is known as the Bethe variational problem (BVP) (Wainwright & Jordan, 2008).
Such a problem can be solved by coordinate descent, where we alternate between updating q̃ to minimize LBethe(θ, q̃) and updating θ to maximize LBethe(θ, q̃). According to Yedidia et al. (2005), updating q̃ to minimize LBethe(θ, q̃) can be exactly achieved by running sum-product loopy belief propagation on pθ, where a fixed point of the belief propagation algorithm yields a local optima of q̃. On the other hand, updating θ to maximize LBethe(θ, q̃) can be easily achieved by gradient ascent. In addition to that, a stationary point (θ∗, q̃∗) of the above BVP is specified by following conditions:
∂LBethe(θ∗, q̃∗) ∂q̃∗ = 0 ∂LBethe(θ∗, q̃∗) ∂θ∗ = 0. (41)
According to Yedidia et al. (2005) and Wainwright & Jordan (2008), the first condition is equivalent to the condition that q̃∗ is specified by a fixed-point of sum-product loopy belief propagation. The second condition states that the moment-matching conditions are satisfied, i.e., qs(ys) = Iy∗s {ys} on each node and qst(ys, yt) = I(y∗s ,y∗t ){(ys, yt)} on each edge.
For our proposed approach SMN, it can be viewed as solving the BVP as defined in Eq. (39). The pre-training stage of SMN aims at initializing θ at a state where the conditions of stationary points in Eq. (41) are approximately satisfied. Then the fine-tuning stage of SMN further adjusts θ to solve the maximin game by alternatively updating θ and q̃.
More specifically, in the pre-training stage, by initializing θ in the way defined by Eq. (27), the collection of pseudomarginal distributions {τs(ys)}s∈V and {τst(ys, yt)}(s,t)∈E is specified by a fixed point of sum-product loopy belief propagation according to proposition 3.1. This implies that ∂∂q̃LBethe(θ, q̃) = 0 for q̃ = {τs(ys)}s∈V ∪ {τst(ys, yt)}(s,t)∈E . Meanwhile, as {τs}s∈V and {τst}(s,t)∈E are learned to match the true labels y∗V on each training graph, we thus have τs(ys) ≈ Iy∗s {ys} on each node and τst(ys, yt) ≈ I(y∗s ,y∗t ){(ys, yt)} on each edge. Therefore, the conditions in Eq. (41) are approximately satisfied by (θ, q̃) with q̃ = {τs(ys)}s∈V ∪ {τst(ys, yt)}(s,t)∈E , which means that the pre-training stage yields a θ to roughly match the conditions of stationary points for the BVP in Eq. (39). After pre-training, the fine-tuning stage of SMN is exactly trying to solve the maximin game of BVP in Eq. (39), where the fine-tuning stage alternates between updating q̃ to minimize LBethe(θ, q̃) via sum-product loopy belief propagation and updating θ to maximize LBethe(θ, q̃) via gradient ascent. As a result, we see that the SMN can also be understood as solving the Bethe variational problem in Eq. (39), which acts as a surrogate for the log-likelihood function.

heading: F EXPERIMENTAL DETAILS
text: Next, we describe our experimental setup in more details.
F.1 DATASET STATISTICS
The statistics of the datasets used in our experiment are summarized in Tab. 8. For the Cora, Citeseer, Pubmed, and PPI datasets, they are under the MIT license.
F.2 MORE DETAILS ABOUT THE DBLP DATASET
Due to the limited space, the experiment section only briefly describes the DBLP dataset, and we complement the details of the dataset here.
The DBLP dataset is constructed from the citation network 5 in Tang et al. (2008). Scientific papers from eight conferences are treated as nodes, and we divide them into three categories for classification based on conference domains 6. For each paper, we compute the mean GloVe embedding 7 (Pennington et al., 2014) of words in the title and abstract as node features. We split the dataset into three disjoint graphs for training/validation/test. The training graph contains papers published before 1999 (with 1999 included). The validation graph contains papers published between 2000 and 2009 (with 2000 and 2009 included). The test graph contains papers published after 2010 (with 2010 included).
5https://originalstatic.aminer.cn/misc/dblp.v12.7z 6ML: ICML/NeurIPS. CV: ICCV/CVPR/ECCV. NLP: ACL/EMNLP/NAACL. 7http://nlp.stanford.edu/data/glove.6B.zip
There exists an undirected edge between two papers if one paper cites the other one. Cross-split edges (e.g., an edge between a paper in the training set and a paper in the validation set) are removed.
F.3 ARCHITECTURE CHOICES
To facilitate reproducibility, we use the GNN module implementations of PyTorch Geometric (Fey & Lenssen, 2019), and follow the GNN models provided in the examples of the repository, unless otherwise mentioned. Note that most architecture choices are not optimal on the benchmark datasets, but we did not tune them since we only need to show that our method brings consistent and significant improvement.
GCN (Kipf & Welling, 2017). We set the number of hidden neurons to 16, and the number of layers to 2. ReLU (Nair & Hinton, 2010) is used as the activation function. We do not dropout between GNN layers.
GraphSage (Hamilton et al., 2017). We set the number of hidden neurons to 64, and the number of layers to 2. ReLU (Nair & Hinton, 2010) is used as the activation function. We do not dropout between GNN layers.
GAT (Veličković et al., 2018). We set the number of hidden neurons to 256 per attention head, and the number of layers to 3. The number of heads for each layer is set to 4, 4 and 6. ELU (Clevert et al., 2016) is used as the activation function. We do not dropout between GNN layers.
Graph U-Net (Gao & Ji, 2019). We set the number of hidden neurons to 64 and the number of layers to 3. We randomly dropout 20% of the edges from the adjacency matrix. We do not dropout node features or between layers.
GCNII (Chen et al., 2020a). We set the number of hidden neurons to 2048 for the citation datasets (Cora*, Citeseer*, Pubmed* and DBLP) and 256 for the PPI dataset. We set the number of layers to 9. ReLU (Nair & Hinton, 2010) is used as the activation function. For PPI, layer normalization (Ba et al., 2016) is applied between the GCNII layers. We do not dropout between GNN layers. We set the strength α of the initial residual connection to 0.5, and the hyperparameter θ to compute the strength of the identity mapping to 1.
The g function. In Eq. (8) of the model section, we define g as a function mapping a pair of Ldimensional representations to a (|Y| × |Y|)-dimensional logit. Two variants of this function are used in our experiment. For the PPI and DBLP dataset, we use the linear variant, where the pair of node representations are concatenated and plugged into a linear layer:
glinear(vs,vt) = W[vs;vt] + b, (42)
where W ∈ R(|Y|×|Y|)×2L is the weight matrix and b ∈ R|Y|×|Y| is the bias. For the citation datasets (Cora*, Citeseer*, Pubmed*), we use the bilienar variant, where the pair of node representations are plugged in a bilinear mapping:
gbilinear(vs,vt) = (Wvs)(Wvt) T , (43)
where W ∈ R|Y|×L is a weight matrix.
SMN with a shared GNN. By default, the SMN uses a node GNN and an edge GNN to approximate the pseudomarginals on nodes and edges respectively. In the experiment, we also consider using a shared GNN for both pseudomarginals on nodes and edges. In other words, us = vs,∀s ∈ V (see Eq. (7) and Eq. (8)). All the other components are the same as the default SMN. The results of this variant are shown in Tab. 6 of the experiment section.
F.4 HYPERPARAMETER CHOICES
GNNs and SMNs. For the node classification task, the learning rate of the node GNN τs in GNNs and SMNs is presented in Tab. 9. For the edge classification task, the learning rate of the edge GNN τst is presented in Tab. 10. For the temperature γ used in the edge GNN τst of SMNs, we summarize its values in Tab. 11.
CRF-linear. For CRF-linear training, we set the learning rate to 5× 10−4.
CRF-GNNs and SMN. For CRF and the refinement stage of SMN, we set the learning rate to 1× 10−5.
GMNN. For GMNN training, we set the learning rate to 5× 10−3.
F.5 COMPUTATIONAL RESOURCES
We run the experiment by using NVIDIA Tesla V100 GPUs with 16GB memory.

heading: G ADDITIONAL RESULTS
text: In this section, we present some additional experimental results.
G.1 ADDITIONAL ANALYSIS OF GNN ARCHITECTURES
In this analysis, we study the effect of node/edge GNN architectures on SMNs. We fix one of the GNNs and change the capacity of the other (Fig. 4), then evaluate SMN-GAT on PPI-1-0, a subset of PPI-1 that only contains its first label. The results show that our model benefit from capacity gain in both node and edge GNNs, highlighting their effective synergy. This also explains the underperformance of SMN-GCN in Tab. 1, where the edge GCN backbone with only two layers and 16 hidden neurons is incapable of modeling the edge label dependencies and thus drags the performance behind. We also find that the node and edge GNNs need not share the same backbone, and in many cases SMNs with different node and edge GNNs perform superior to those with same backbone (Fig. 4). The expressiveness of edge GNNs is crucial to the performance of SMN. Though we did not optimize the design of our edge GNNs, they have shown to be helpful in boosting the performance once plugged into our approach.
G.2 NODE-LEVEL ACCURACY ON CORA*, CITESEER*, AND PUBMED*
In the experiment, we report the graph-level accuracy on the Cora*, Citeseer*, and Pubmed* datasets, where SMNs consistently outperform other methods. Besides the graph-level accuracy, we also
compute the node-level accuracy on these datasets, and the results are reported in Tab. 12. We can see that our approach still consistently outperforms other methods in terms of node-level accuracy.
G.3 COMPARISON OF SUM-PRODUCT AND MAX-PRODUCT BELIEF PROPAGATION
As explained in section 4.2, the sum-product belief propagation algorithm is more applicable to the case of node-level accuracy, as it aims at inferring the marginal label distribution on each node. Nevertheless, in practice we find that the max-product algorithm usually achieves better empirical node-level accuracy. For example, the results on the PPI-10 dataset are presented in Tab. 13.
Because of the better empirical results, we choose to use max-product belief propagation by default.
G.4 HYPERPARAMETER ANALYSIS
Finally, We present analysis of the hyperparameter γ (i.e., edge temperature) in Fig. 5.

Review summaries:
The paper targets the task of graph node classification in the inductive setting, taking an input node and edge feature representations, and inferring a categorical label for each node. It focuses on the case where it is not sufficient to classify a node independently of its neighbours, but where information stemming from predictions for neighbouring nodes needs to be taken into account.  It improves over the current approaches by offering a performant and efficient method to combine advantages of CRFs (joint inference) with the representational power afforded by GNNs, while overcoming the computational challenges. In particular, this paper improves on the closest attempts at combining CRF and GNN, Ma+ 2018 and Qu+ 2019, which use pseudolikelihood; Qu+ 2019 is used as baseline in experiments (“GMNN”), while pseudoloikelihood as training objective is investigated in sec5.5.5.  --- update after rebuttals I have read other reviews and the authors' replies, as well as extra experiments. I am maintaining my score.

The paper aims to break the label independence assumption in existing GNNs. To this end, a combination of GNNs and CRF is proposed, named SMN. SMN improves GNNs by modeling the joint distribution among discrete node labels; SMN is a promising alternative to CRF as it provides efficient learning and inference procedures. Empirically, SMN achieves better node-level and graph-level prediction accuracies than several existing models, with minor additional runtime overhead.

This paper proposes a CRF for classifying nodes in graphs where the CRF has a potential for each node and edge in the graph. GNNs are used for computing the potentials, one GNN for node potentials and one GNN for edge potentials. For general graphs, computing the partition function is intractable, so approximations are used during both learning and inference.  Learning draws from prior work and combines learning pseudomarginals of nodes and edges with GNNs and optionally some steps of "refinement" by optimizing a maximin game equivalent to likelihood maximization. Inference uses sum-product or max-product loopy BP (they perform similarly, though max-product is slightly better). The procedure is called a "Structured Markov Network" (SMN). Experimental results on several graph node classification tasks show that SMNs outperform several baselines, including both CRFs with standard training (using approximate inference during training) and CRFs trained with pseudolikelihood.  

Authors study the problem of node labeling in the inductive case, i.e., at test time the goal is to label all the nodes of a given graph. For that problem, several variants of GNNs and CRFs have been proposed in the past, and the authors propose a Structural Markov Network that uses GNNs to model the potential functions of a CRF with the subtle difference that a proxy optimization problem is solved to make learning more efficient. Experiments are provided to demonstrate the applicability of their method.