A Quasi-Wasserstein Loss for Learning Graph Neural Networks
When learning graph neural networks (GNNs) in node-level prediction tasks,
most existing loss functions are applied for each node independently, even if
node embeddings and their labels are non-i.i.d. because of their graph
structures. To eliminate such inconsistency, in this study we propose a novel
Quasi-Wasserstein (QW) loss with the help of the optimal transport defined on
graphs, leading to new learning and prediction paradigms of GNNs. In
particular, we design a ``Quasi-Wasserstein'' distance between the observed
multi-dimensional node labels and their estimations, optimizing the label
transport defined on graph edges. The estimations are parameterized by a GNN in
which the optimal label transport may determine the graph edge weights
optionally. By reformulating the strict constraint of the label transport to a
Bregman divergence-based regularizer, we obtain the proposed Quasi-Wasserstein
loss associated with two efficient solvers learning the GNN together with
optimal label transport. When predicting node labels, our model combines the
output of the GNN with the residual component provided by the optimal label
transport, leading to a new transductive prediction paradigm. Experiments show
that the proposed QW loss applies to various GNNs and helps to improve their
performance in node-level classification and regression tasks. The code of this
work can be found at \url{https://github.com/SDS-Lab/QW_Loss}.
Minjie Cheng, Hongteng Xu
# A Quasi-Wasserstein Loss for Learning Graph Neural Networks

###### Abstract

When learning graph neural networks (GNNs) in node-level prediction tasks, most existing loss functions are applied for each node independently, even if node embeddings and their labels are non-i.i.d. because of their graph structures. To eliminate such inconsistency, in this study we propose a novel Quasi-Wasserstein (QW) loss with the help of the optimal transport defined on graphs, leading to new learning and prediction paradigms of GNNs. In particular, we design a "Quasi-Wasserstein" distance between the observed multi-dimensional node labels and their estimations, optimizing the label transport defined on graph edges. The estimations are parameterized by a GNN in which the optimal label transport may determine the graph edge weights optionally. By reformulating the strict constraint of the label transport to a Bregman divergence-based regularizer, we obtain the proposed Quasi-Wasserstein loss associated with two efficient solvers learning the GNN together with optimal label transport. When predicting node labels, our model combines the output of the GNN with the residual component provided by the optimal label transport, leading to a new transductive prediction paradigm. Experiments show that the proposed QW loss applies to various GNNs and helps to improve their performance in node-level classification and regression tasks.

## 1 Introduction

Graph neural network (GNN) plays a central role in many graph learning tasks, such as social network analysis (Qiu et al., 2018; Fan et al., 2019; Zhang et al., 2023), molecular modeling (Satorras et al., 2021; Jiang et al., 2021; Wang et al., 2022), transportation forecasting (Wang et al., 2020; Li and Zhu, 2021), and so on. Given a graph with node features, a GNN embeds the graph nodes by exchanging and aggregating the node features, whose implementation is based on message-passing operators in the spatial domain (Niepert et al., 2016; Kipf and Welling, 2017; Velickovic et al., 2018) or graph filters in the spectral domain (Defferrard et al., 2016; Chien et al., 2021; Bianchi et al., 2021; He et al., 2021). When some node labels are available, we can learn the GNN in a node-level semi-supervised learning (Yang et al., 2016; Kipf and Welling, 2017; Xu et al., 2019), optimizing the node embeddings to predict the observed labels. This learning framework has achieved encouraging performance in many node-level prediction tasks, e.g., node classification (Sen et al., 2008; McAuley et al., 2015). When applying the above node-level GNN learning framework, existing work often leverages a loss function (e.g., the cross-entropy loss) to penalize the discrepancy between each
node's label and the corresponding estimation. Here, some inconsistency between the objective design and the intrinsic data structure arises -- the objective of learning a GNN is implemented as the summation of all the nodes' loss functions, which is often applied for i.i.d. data, but the node embeddings and their labels are non-i.i.d. in general because of the underlying graph structure and the information aggregation achieved by the GNN. As a result, the current objective treats the losses of individual nodes independently and evenly, even if the nodes in the graph are correlated and have different significance for learning the GNN. Such inconsistency may lead to sub-optimal GNNs in practice, but to our knowledge none of existing work considers this issue in-depth.

To eliminate the inconsistency, we leverage computational optimal transport techniques (Peyre and Cuturi 2019), proposing a new objective called Quasi-Wasserstein (QW) loss for learning GNNs. As illustrated in Figure 1, given partially-observed node labels and their estimations parametrized by a GNN, we consider the optimal transport between them and formulate the problem as the aggregation of the Wasserstein distances (Frogner et al. 2015) corresponding to all label dimensions. This problem can be equivalently formulated as a label transport minimization problem (Essid and Solomon 2018, Facca and Benzi 2021) defined on the graph, leading to the proposed QW loss. By minimizing this loss, we can jointly learn the optimal label transport and the GNN parametrizing the label estimations. This optimization problem can be solved efficiently by Bregman divergence-based algorithms, e.g., Bregman ADMM (Wang and Banerjee 2014, Xu 2020). Optionally, through

Figure 1: The scheme of our QW-loss and the corresponding learning paradigm. Given a graph whose node features are denoted as blue circles and partially-observed node labels are denoted as blue stems, a GNN embeds the graph nodes and outputs estimated labels (denoted as orange stems). By minimizing the QW loss, we obtain the optimal label transport (denoted as the dotted red arrows on the graph edges) between the real and estimated node labels. Optionally, the optimal label transport can be used to determine the weights of graph edges (through an edge weight predictor). The final predictions are the combinations of the optimal label transport and the estimated node labels.
a multi-layer perceptron (MLP), we can determine the edge weights of the graph based on optimal label transport, leading to a GNN with learnable edge weights.

The contributions of this study include the following two points:

* **A theoretically-solid loss without the inconsistency issue.** The QW loss provides a new optimal transport-based loss for learning GNNs, which considers the labels and estimations of graph nodes jointly. Without the questionable i.i.d. assumption, it eliminates the inconsistency issue mentioned above. In theory, we demonstrate that the QW loss is a valid metric for the node labels defined on graphs. Additionally, the traditional objective function for learning GNNs can be treated as a special case of our QW loss. We further demonstrate that applying our QW loss reduces data fitting errors in the training phase.
* **New learning and prediction paradigms.** Different from the existing methods that combine GNNs with label propagation mechanisms (Huang et al., 2021; Wang and Leskovec, 2021; Dong et al., 2021), the QW loss provides a new way to combine node embeddings with label information in both training and testing phases. In particular, Bregman divergence-based algorithms are applied to learn the model, and the final model consists of the GNN and a residual component provided by the optimal label transport. When predicting node labels, the model combines the estimations provided by the GNN with the complementary information from the optimal label transport, leading to a new transductive prediction paradigm.

Experiments demonstrate that our QW loss applies to various GNNs and helps to improve their performance in various node-level classification and regression tasks.

## 2 Related work

### Graph Neural Networks

Graph neural networks can be coarsely categorized into two classes. The GNNs in the first class apply spatial convolutions to graphs (Niepert et al., 2016). The representative work includes the graph convolutional network (GCN) in (Kipf and Welling, 2017), the graph attention network (GAT) in (Velickovic et al., 2018), and their variants (Zhuang and Ma, 2018; Wang et al., 2019; Fu et al., 2020). The GNNs in the second class achieve graph spectral filtering (Guo et al., 2023). They are often designed based on a polynomial basis, such as ChebNet (Defferrard et al., 2016) and its variants (He et al., 2022), GPR-GNN (Chien et al., 2021), and BernNet (He et al., 2021). Besides approximated by the polynomial basis, the spectral GNNs can be learned by other strategies, e.g., Personalized PageRank in APPNP (Gasteiger et al., 2019), graph optimization functions in GNN-LF/HF (Zhu et al., 2021), ARMA filters (Bianchi et al., 2021), and diffusion kernel-based filters (Klicpera et al., 2019; Xu et al., 2019; Du et al., 2022).

The above spatial and spectral GNNs are correlated because a spatial convolution always corresponds to a graph spectral filter (Balcilar et al., 2021). For example, GCN (Kipf and Welling, 2017) can be explained as a low-pass filter achieved by a first-order Chebyshev polynomial. Given a graph with some labeled nodes, we often learn the above GNNs in a semi-supervised node-level learning framework (Kipf and Welling, 2017; Yang et al., 2016), in which the GNNs embed all the nodes and are trained under the supervision of the labeled nodes. However, as aforementioned, the objective functions used in the framework treat the graph nodes independently and thus mismatch with the non-i.i.d. nature of the data.
### Computational Optimal Transport

As a powerful mathematical tool, optimal transport (OT) distance (or called Wasserstein distance under some specific settings) provides a valid metric for probability measures (Villani 2008), which has been widely used for various machine learning problems, e.g., distribution comparison (Frogner et al. 2015, Lee et al. 2019), point cloud registration (Grave et al. 2019), graph partitioning (Dong and Sawin 2020, Xu et al. 2019b), generative modeling (Arjovsky et al. 2017, Tolstikhin et al. 2018), and so on. Typically, the OT distance corresponds to a constrained linear programming problem. To approximate the OT distance with low complexity, many algorithms have been proposed, e.g., Sinkhorn-scaling (Cuturi 2013), Bregman ADMM (Wang and Banerjee 2014), Conditional Gradient (Titouan et al. 2019), and Inexact Proximal Point (Xie et al. 2020). Recently, two iterative optimization methods have been proposed to solve the optimal transport problems defined on graphs (Essid and Solomon 2018, Facca and Benzi 2021).

These efficient algorithms make the OT distance a feasible loss for machine learning problems, e.g., the Wasserstein loss in (Frogner et al. 2015). Focusing on the learning of GNNs, the work in (Chen et al. 2020) proposes a Wasserstein distance-based contrastive learning method. The Gromovized Wasserstein loss is applied to learn cross-graph node embeddings (Xu et al. 2019b), graph factorization models (Vincent-Cuaz et al. 2021, Xu 2020), and GNN-based graph autoencoders (Xu et al. 2022). The above work is designed for graph-level learning tasks, e.g., graph matching, representation, classification, and clustering. Our QW loss, on the contrary, is designed for node-level prediction tasks, resulting in significantly different learning and prediction paradigms.

## 3 Proposed Method

### Motivation and Principle

Denote a graph as \(G(\mathcal{V},\mathcal{E})\), where \(\mathcal{V}\) represents the set of nodes and \(\mathcal{E}\) represents the set of edges, respectively. The graph \(G\) is associated with an adjacency matrix \(\mathbf{A}\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}\) and a edge weight vector \(\mathbf{w}=[w_{e}]\in\mathbb{R}^{|\mathcal{E}|}\). The weights in \(\mathbf{w}\) correspond to the non-zero elements in \(\mathbf{A}\). For an unweighted graph, \(\mathbf{A}\) is a binary matrix, and \(\mathbf{w}\) is an all-one vector. Additionally, the nodes of the graph may have \(D\)-dimensional features, which are formulated as a matrix \(\mathbf{X}\in\mathbb{R}^{|\mathcal{V}|\times D}\). Suppose that a subset of nodes, denoted as \(\mathcal{V}_{L}\subset\mathcal{V}\), are annotated with \(C\)-dimensional labels, i.e., \(\{\mathbf{y}_{v}\in\mathbb{R}^{C}\}_{v\in\mathcal{V}_{L}}\). We would like to learn a GNN to predict the labels of the remaining nodes, i.e., \(\{\mathbf{y}_{v}\}_{v\in\mathcal{V}\setminus\mathcal{V}_{L}}\).

The motivation for applying GNNs is based on the non-i.i.d. property of the node features and labels. Suppose that we have two nodes connected by an edge, i.e., \((v,v^{\prime})\in\mathcal{E}\), where \((\mathbf{x}_{v},\mathbf{y}_{v})\) and \((\mathbf{x}_{v^{\prime}},\mathbf{y}_{v^{\prime}})\) are their node features and labels. For each node, its neighbors' features or labels can provide valuable information to its prediction task, i.e., the conditional probability \(p(\mathbf{y}_{v}|\mathbf{x}_{v})\neq p(\mathbf{y}_{v}|\mathbf{x}_{v},\mathbf{x}_{v^{\prime}})\) and \(p(\mathbf{y}_{v}|\mathbf{x}_{v})\neq p(\mathbf{y}_{v}|\mathbf{x}_{v},\mathbf{y}_{v^{\prime}})\) in general. Similarly, for node pairs, their labels are often conditionally-dependent, i.e., \(p(\mathbf{y}_{v},\mathbf{y}_{v^{\prime}}|\mathbf{x}_{v},\mathbf{x}_{v^{\prime}})=p(\mathbf{y}_{v}| \mathbf{x}_{v},\mathbf{x}_{v^{\prime}},\mathbf{y}_{v^{\prime}})p(\mathbf{y}_{v^{\prime}}|\mathbf{ x}_{v},\mathbf{x}_{v^{\prime}})\neq p(\mathbf{y}_{v}|\mathbf{x}_{v},\mathbf{x}_{v^{\prime}})p(\mathbf{y}_{v ^{\prime}}|\mathbf{x}_{v},\mathbf{x}_{v^{\prime}})\). More generally, for all node labels, we have

\[p(\{\mathbf{y}_{v}\}_{v\in\mathcal{V}}|\mathbf{X},\mathbf{A})\neq\,\prod\nolimits_{v\in \mathcal{V}}p(\mathbf{y}_{v}|\mathbf{X},\mathbf{A}), \tag{1}\]

Ideally, we shall learn a GNN to maximize the conditional probability of all labeled nodes, i.e., \(\max p(\{\mathbf{y}_{v}\}_{v\in\mathcal{V}_{L}}|\mathbf{X},\mathbf{A})\). In practice, however, most existing methods formulate the node-level learning paradigm of the GNN as

\[\max_{\theta}\prod\nolimits_{v\in\mathcal{V}_{L}}p(\mathbf{y}_{v}|\mathbf{X},\mathbf{A} ;\theta)\Leftrightarrow\,\min_{\theta}\,\sum\nolimits_{v\in\mathcal{V}_{L}} \psi(g_{v}(\mathbf{X},\mathbf{A};\theta),\ \mathbf{y}_{v}). \tag{2}\]
Here, \(g\) is a graph neural network whose parameters are denoted as \(\theta\). Taking the adjacency matrix \(\mathbf{A}\) and the node feature matrix \(\mathbf{X}\) as input, the GNN \(g\) predicts the node labels. \(g_{v}(\mathbf{X},\mathbf{A};\theta)\) represents the estimation of the node \(v\)'s label achieved by the GNN, which is also denoted as \(\hat{\mathbf{y}}_{v}\). Similarly, we denote \(g_{\mathcal{V}}(\mathbf{X},\mathbf{A};\theta)\) as the estimated labels for the node set \(\mathcal{V}\) in the following content. The loss function \(\psi:\mathbb{R}^{C}\times\mathbb{R}^{C}\mapsto\mathbb{R}\) is defined in the node level. In node-level classification tasks, it is often implemented as the cross-entropy loss or the KL-divergence (i.e., \(p(\mathbf{y}_{v}|\mathbf{X},\mathbf{A};\theta)\) is modeled by the softmax function). In node-level regression tasks, it is often implemented as the least-square loss (i.e., \(p(\mathbf{y}_{v}|\mathbf{X},\mathbf{A};\theta)\) is assumed to be the Gaussian distribution).

**The loss in (2) assumes the node labels to be conditionally-independent with each other, which may be too strong in practice and inconsistent with the non-i.i.d. property of graph-structured data shown in (1).** To eliminate such inconsistency, we should treat node labels as a set rather than independent individuals, developing a set-level loss to penalize the discrepancy between the observed labels and their estimations globally, i.e.,

\[\min_{\theta}\text{Loss}(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta),\{\mathbf{y}_{ v}\}_{v\in\mathcal{V}_{L}}). \tag{3}\]

In the following content, we will design such a loss with theoretical supports, based on the optimal transport on graphs.

### Optimal Transport on Graphs

Suppose that we have two measures on a graph \(G(\mathcal{V},\mathcal{E})\), denoted as \(\mathbf{\mu}\in[0,\infty)^{|\mathcal{V}|}\) and \(\mathbf{\gamma}\in[0,\infty)^{|\mathcal{V}|}\), respectively. The element of each measure indicates the "mass" of a node. Assume the two measures to be balanced, i.e., \(\langle\mathbf{\gamma}-\mathbf{\mu},\mathbf{1}_{|\mathcal{V}|}\rangle=0\), where \(\langle\cdot,\cdot\rangle\) is the inner product operator. The optimal transport, or called the 1-Wasserstein distance (Villani, 2008), between them is defined as

\[W_{1}(\mathbf{\mu},\mathbf{\gamma}):=\,\min_{\mathbf{T}\in\Pi(\mathbf{\mu},\mathbf{\gamma})} \langle\mathbf{D},\;\mathbf{T}\rangle=\,\min_{\mathbf{T}\in\Pi(\mathbf{\mu},\mathbf{\gamma})}\, \sum\nolimits_{v,v^{\prime}\in\mathcal{V}\times\mathcal{V}}t_{vv^{\prime}}d_ {vv^{\prime}}, \tag{4}\]

where \(\mathbf{D}=[d_{vv^{\prime}}]\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}\) represents the shortest path distance matrix, and \(\Pi(\mathbf{\mu},\mathbf{\gamma})=\{\mathbf{T}\geq\mathbf{0}|\mathbf{T}\mathbf{1}_{|\mathcal{V}|}=\bm {\mu},\mathbf{T}^{\top}\mathbf{1}_{|\mathcal{V}|}=\mathbf{\gamma}\}\) represents the set of all valid doubly stochastic matrices. Each \(\mathbf{T}=[t_{vv^{\prime}}]\in\Pi(\mathbf{\mu},\mathbf{\gamma})\) is a transport plan matrix. The optimization problem in (4) corresponds to finding the optimal transport plan \(\mathbf{T}^{*}=[t_{vv^{\prime}}^{*}]\) to minimize the "cost" of changing \(\mathbf{\mu}\) to \(\mathbf{\gamma}\), in which the cost is measured as the sum of "mass" \(t_{vv^{\prime}}\) moved from node \(v\) to node \(v^{\prime}\) times distance \(d_{vv^{\prime}}\).

#### 3.2.1 Wasserstein Distance for Vectors on A Graph

For the optimal transport problem defined on graphs, we can simplify the problem in (4) by leveraging the underlying graph structures. As shown in (Essid and Solomon, 2018; Santambrogio, 2015), given a graph \(G(\mathcal{V},\mathcal{E})\), we can define a sparse matrix \(\mathbf{S}_{\mathcal{V}}=[s_{ve}]\in\{0,\pm 1\}^{|\mathcal{V}|\times|\mathcal{E}|}\) to indicating the graph topology. For node \(v\) and edge \(e\), the corresponding element in \(\mathbf{S}_{\mathcal{V}}\) is

\[s_{ve}=\begin{cases}1&\text{if $v$ is ``head'' of edge $e$}\\ -1&\text{if $v$ is ``tail'' of edge $e$}\\ 0&\text{otherwise.}\end{cases} \tag{5}\]

When \(G\) is directed, the "head" and "tail" of each edge are predefined. When \(G\) is undirected, we can randomly define each edge's "head" and "tail". Accordingly, the 1-Wasserstein distance in (4) can be equivalently formulated as a minimum-cost flow problem:

\[W_{1}(\mathbf{\mu},\mathbf{\gamma})=\,\min_{\mathbf{f}\in\Omega(\mathbf{S}_{\mathcal{V}},\mathbf{ \mu},\mathbf{\gamma})}\|\text{diag}(\mathbf{w})\mathbf{f}\|_{1}, \tag{6}\]
where \(\mathrm{diag}(\mathbf{w})\) is a diagonal matrix constructed by the edge weights of the graph. The vector \(\mathbf{f}=[f_{e}]\in\Omega(\mathbf{S}_{\mathcal{V}},\mathbf{\mu},\mathbf{\gamma})\) is the flow indicating the mass passing through each edge. Accordingly, the cost corresponding to edge \(e\) is represented as the distance \(w_{e}\) times the mass \(f_{e}\), and the sum of all the costs leads to the objective function in (6). The feasible domain of the flow vector is defined as

\[\Omega(\mathbf{S}_{\mathcal{V}},\mathbf{\mu},\mathbf{\gamma})=\mathcal{U}^{| \mathcal{E}|}\cap\{\mathbf{f}\ |\ \mathbf{S}_{\mathcal{V}}\mathbf{f}=\mathbf{\gamma}-\mathbf{\mu}\},\quad\text{ where }\mathcal{U}=\begin{cases}[0,\infty),&\text{Directed }G,\\ \mathbb{R},&\text{Undirected }G.\end{cases} \tag{7}\]

The constraint \(\mathbf{S}_{\mathcal{V}}\mathbf{f}=\mathbf{\gamma}-\mathbf{\mu}\) ensures that the flow on all the edges leads to the change from \(\mathbf{\mu}\) to \(\mathbf{\gamma}\). The flow in the directed graph is nonnegative, which can only pass through each edge from "head" to "tail". On the contrary, the undirected graph allows the mass to transport from "tail" to "head". By solving (6), we find the optimal flow \(\mathbf{f}^{*}\) (or the so-called optimal transport on edges) that minimizes the overall cost \(\|\mathrm{diag}(\mathbf{w})\mathbf{f}\|_{1}\).

#### 3.2.2 Partial Wasserstein Distance for Vectors on a Graph

When only a subset of nodes (i.e., \(\mathcal{V}_{L}\subset\mathcal{V}\)) have observable signals, we can define a _partial Wasserstein distance_ on the graph by preserving the constraints relevant to \(\mathcal{V}_{L}\): for \(\mathbf{\mu},\mathbf{\gamma}\in\mathbb{R}^{|\mathcal{V}_{L}|}\), we have

\[W_{1}^{(P)}(\mathbf{\mu},\mathbf{\gamma})=\,\min\nolimits_{\mathbf{f}\in \Omega(\mathbf{S}_{\mathcal{V}_{L}},\mathbf{\mu},\mathbf{\gamma})}\|\mathrm{diag}(\mathbf{w}) \mathbf{f}\|_{1}, \tag{8}\]

where \(\mathbf{S}_{\mathcal{V}_{L}}\in\{0,\pm 1\}^{|\mathcal{V}_{L}|\times|\mathcal{E}|}\) is a submatrix of \(\mathbf{S}_{\mathcal{V}}\), storing the rows corresponding to the nodes with observable signals.

Different from the classic Wasserstein distance, the Wasserstein distance defined on graphs is not limited for nonnegative vectors because the minimum-cost flow formulation in (6) is _shift-invariant_, i.e., \(\forall\mathbf{\delta}\in\mathbb{R}^{|\mathcal{V}|}\), \(W_{1}(\mathbf{\mu},\mathbf{\gamma})=W_{1}(\mathbf{\mu}-\mathbf{\delta},\mathbf{\gamma}-\mathbf{\delta})\). The partial Wasserstein distance in (8) can be viewed as a generalization of (6), and thus it holds the shift-invariance as well. Moreover, given an undirected graph, we can prove that both the \(W_{1}\) in (6) and the \(W_{1}^{(P)}\) in (8) can be valid metrics in the spaces determined by the topology of an undirected graph (See Appendix A for a complete proof).

**Theorem 1** (The Validness as A Metric).: _Given an undirected graph \(G(\mathcal{V},\mathcal{E})\), with edge weights \(\mathbf{w}\in[0,\infty)^{|\mathcal{E}|}\) and a matrix \(\mathbf{S}_{\mathcal{V}}\in\{0,\pm 1\}^{|\mathcal{V}|\times\mathcal{E}}\) defined in (5), the \(W_{1}\) in (6) is a metric in \(Range(\mathbf{S}_{\mathcal{V}})\), and the \(W_{1}^{(P)}\) in (8) is a metric in \(Range(\mathbf{S}_{\mathcal{V}_{L}})\), \(\forall\mathcal{V}_{L}\in\mathcal{V}\)._

The flow vector \(\mathbf{f}\) in (8) has fewer constraints than that in (6), so the relation between \(W_{1}\) and \(W_{1}^{(P)}\) obeys the following theorem.

**Theorem 2** (The Monotonicity).: _Given an undirected graph \(G(\mathcal{V},\mathcal{E})\), with edge weights \(\mathbf{w}\in[0,\infty)^{|\mathcal{E}|}\) and a matrix \(\mathbf{S}_{\mathcal{V}}\in\{0,\pm 1\}^{|\mathcal{V}|\times\mathcal{E}}\) defined in (5), we have_

\[W_{1}(\mathbf{\mu},\mathbf{\gamma})\geq W_{1}^{(P)}(\mathbf{\mu}_{\mathcal{V}^{\prime}}, \mathbf{\gamma}_{\mathcal{V}^{\prime\prime}})\geq W_{1}^{(P)}(\mathbf{\mu}_{\mathcal{ V}^{\prime\prime}},\mathbf{\gamma}_{\mathcal{V}^{\prime\prime}})\geq 0,\quad\forall\mathbf{\mu},\mathbf{ \gamma}\in Range(\mathbf{S}_{\mathcal{V}}),\ \mathcal{V}^{\prime\prime}\subset\mathcal{V} ^{\prime}\subset\mathcal{V}, \tag{9}\]

_where \(\mathbf{\mu}_{\mathcal{V}^{\prime}}\) denotes the subvector of \(\mathbf{\mu}\) corresponding to the set \(\mathcal{V}^{\prime}\)._

### Learning GNNs with Quasi-Wasserstein Loss

Inspired by the optimal transport on graphs, we propose a Quasi-Wasserstein loss for learning GNNs. Given a partially-labeled graph \(G\), we formulate observed labels as a matrix
\(\mathbb{R}^{|\mathcal{V}_{L}|\times C}\), where \(\mathbf{y}^{(c)}_{\mathcal{V}_{L}}\in\mathbb{R}^{|\mathcal{V}_{L}|}\) represents the labels in \(c\)-th dimension. The estimated labels generated by a GNN \(g\) are represented as \(\widehat{\mathbf{Y}}_{\mathcal{V}_{L}}=[\hat{\mathbf{y}}^{(c)}_{\mathcal{V}_{L}}]:=g_{ \mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta)\). Our QW loss is defined as

\[\begin{split} QW(\widehat{\mathbf{Y}}_{\mathcal{V}_{L}},\mathbf{Y}_{ \mathcal{V}_{L}})&=\,\sum\nolimits_{c=1}^{C}W_{1}^{(P)}( \hat{\mathbf{y}}^{(c)}_{\mathcal{V}_{L}},\mathbf{y}^{(c)}_{\mathcal{V}_{L}})\\ &=\,\sum\nolimits_{c=1}^{C}\,\min_{\mathbf{f}^{(c)}\in\Omega(\mathbf{S} _{\mathcal{V}_{L}},\,\,\hat{\mathbf{y}}^{(c)}_{\mathcal{V}_{L}},\,\,\hat{\mathbf{y}}^{ (c)}_{\mathcal{V}_{L}})}\|\text{diag}(\mathbf{w})\mathbf{f}^{(c)}\|_{1}\\ &=\,\min_{\mathbf{F}\in\Omega_{C}(\mathbf{S}_{\mathcal{V}_{L}},\,\,g_{ \mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta),\,\,\mathcal{Y}_{\mathcal{V}_{L}})}\| \text{diag}(\mathbf{w})\mathbf{F}\|_{1},\end{split} \tag{10}\]

where \(\mathbf{F}=[\mathbf{f}^{(c)}]\) is the flow matrix, and its feasible domain \(\Omega_{C}=\mathcal{U}^{|\mathcal{E}|\times C}\cap\{\mathbf{F}|\mathbf{S}_{\mathcal{ V}_{L}}\mathbf{F}=\mathbf{Y}_{\mathcal{V}_{L}}-g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A}; \theta)\}\). \(W_{1}^{(P)}(\hat{\mathbf{y}}^{(c)}_{\mathcal{V}_{L}},\mathbf{y}^{(c)}_{\mathcal{V}_{L}})\) is the partial Wasserstein distance between the observed and estimated labels in the \(c\)-th dimension. The QW loss is the summation of the \(C\) partial Wasserstein distances, and it can be rewritten as the summation of \(C\) minimum-cost flow problems. As shown in the third row of (10), these \(C\) problems can be formulated as a single optimization problem, and the variables of the problems are aggregated as the flow matrix \(\mathbf{F}\). Based on Theorem 1, our QW loss is a metric for the matrices whose columns are in \(Range(\mathbf{S}_{\mathcal{V}_{L}})\).

The optimal flow matrix \(\mathbf{F}^{*}=[\mathbf{f}^{(c)*}]\) is called **optimal label transport** in this study, in which the column \(\mathbf{f}^{(c)*}\) indicating the optimal transport on graph edges between the observed labels and their estimations in the \(c\)-th dimension. Note that, we call the proposed loss "Quasi-Wasserstein" because \(i)\) it is not equal to the classic 1-Wasserstein distance between the label sets \(\widehat{\mathbf{Y}}_{\mathcal{V}_{L}}\) and \(\mathbf{Y}_{\mathcal{V}_{L}}\)(Peyre and Cuturi, 2019), and \(ii)\) the optimal label transport \(\mathbf{F}^{*}\) cannot be reformulated as the optimal transport \(\mathbf{T}^{*}\) obtained by (4). When \(\mathbf{F}^{*}=\mathbf{0}\), we have \(QW(\widehat{\mathbf{Y}}_{\mathcal{V}_{L}},\mathbf{Y}_{\mathcal{V}_{L}})=0\) and accordingly, \(\widehat{\mathbf{Y}}_{\mathcal{V}_{L}}=\mathbf{Y}_{\mathcal{V}_{L}}\), which means that the GNN \(g\) perfectly estimates the observed labels. Therefore, the QW loss provides a new alternative for the training loss of GNN. Different from the traditional loss function in (2), the QW loss treats observed node labels and their estimations as two sets and measures their discrepancy accordingly, which provides an effective implementation of the loss in (3).

Applying the QW loss to learn a GNN \(g\) results in the following constrained optimization problem:

\[\min_{\theta}QW(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta),\mathbf{Y}_{\mathcal{V}_ {L}})=\,\min_{\theta}\,\min_{\mathbf{F}\in\Omega_{C}(\mathbf{S}_{\mathcal{V}_{L}},\,g_ {\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta),\,\,\mathbf{Y}_{\mathcal{V}_{L}})}\|\text{ diag}(\mathbf{w})\mathbf{F}\|_{1}. \tag{11}\]

To solve it effectively, we consider the following two algorithms.

#### 3.3.1 Bregman Divergence-based Approximate Solver

By relaxing the constraint \(\{\mathbf{F}|\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}=\mathbf{Y}_{\mathcal{V}_{L}}-g_{\mathcal{V} _{L}}(\mathbf{X},\mathbf{A};\theta)\}\) to a Bregman divergence-based regularizer (Benamou et al., 2015), we can reformulate (11) to the following problem:

\[\begin{split}\min_{\theta,\,\,\mathbf{F}\in\mathcal{U}^{|\mathcal{E}| \times C}}\|\text{diag}(\mathbf{w})\mathbf{F}\|_{1}+\lambda\underbrace{B_{\phi}(g_{ \mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta)+\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F},\,\, \mathbf{Y}_{\mathcal{V}_{L}})}_{\sum_{v\in\mathcal{V}_{L}}\psi(g_{v}(\mathbf{X},\mathbf{A}; \theta)+\mathbf{S}_{v}\mathbf{F},\,\,\mathbf{y}_{v})}.\end{split} \tag{12}\]

Here, \(B_{\phi}(x,y)=\phi(x)-\phi(y)-\langle\nabla\phi(y),x-y\rangle\) is the Bregman divergence defined based on the strictly convex function \(\phi\), and the hyperparameter \(\lambda>0\) controls its significance. In node classification tasks, we can set \(\phi\) as an entropic function, and the Bregman divergence becomes the KL-divergence. In node regression tasks, we can set \(\phi\) as a least-square loss, and the Bregman divergence becomes the least-square loss. Therefore, as shown in (12), the Bregman divergence \(B_{\phi}\) can be implemented based on commonly-used loss functions, e.g., the \(\psi\) in (2).
Algorithm 1 shows the algorithmic scheme in details. For undirected graphs, the \(\mathcal{U}\) in (12) is \(\mathbb{R}\), and accordingly, (12) becomes a unconstrained optimization problem. We can solve it efficiently by gradient descent. For directed graphs, we just need to add a projection step when updating the flow matrix \(\mathbf{F}\), leading to the projected gradient descent algorithm.

#### 3.3.2 Bregman ADMM-based Exact Solver

Algorithm 1 solves (11) approximately -- because of relaxing the strict equality constraint to a regularizer, the solution of (12) often cannot satisfy the original equality constraint. To solve (11) exactly, we further develop a learning method based on the Bregman ADMM (Bremgan Alternating Direction Method of Multipliers) algorithm (Wang and Banerjee 2014). In particular, we can rewrite (11) in the following augmented Lagrangian form:

\[\begin{split}\min_{\theta,\ \mathbf{F}\in\mathcal{U}^{|\mathcal{E}| \times C},\ \mathbf{Z}\in\mathbb{R}^{|\mathcal{V}_{L}|\times C}}&\|\text{diag}( \mathbf{w})\mathbf{F}\|_{1}+\langle\mathbf{Z},\ g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta)+ \mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}-\mathbf{Y}_{\mathcal{V}_{L}}\rangle\\ &+\lambda B_{\phi}(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta)+\mathbf{ S}_{\mathcal{V}_{L}}\mathbf{F},\ \mathbf{Y}_{\mathcal{V}_{L}}),\end{split} \tag{13}\]

where \(\mathbf{Z}\in\mathbb{R}^{|\mathcal{V}_{L}|\times C}\) is the dual variable, the second term in (13) is the Lagrangian term corresponding to the equality constraint, and the third term in (13) is the augmented term implemented as the Bregman divergence.

```
0: A graph \(G\), its adjacency matrix \(\mathbf{A}\) (edge weights \(\mathbf{w}\)), node features \(\mathbf{X}\), observed labels \(\mathbf{Y}_{\mathcal{V}_{L}}\), the number of inner iterations \(J\).
1: Initialize \(\theta^{(0)}\) and \(\mathbf{F}^{(0)}\) randomly and set \(\mathbf{Z}^{(0)}=\mathbf{0}_{|\mathcal{V}_{L}|\times C}\).
2:while Not converge do
3: Solve (14) by Adam (Kingma and Ba 2014) with \(J\) steps and obtain \(\theta^{(k+1)}\).
4: Solve (15) by Adam (Kingma and Ba 2014) with \(J\) steps and obtain \(\mathbf{F}^{(k+1)}\).
5:if\(G\) is a directed graph then\(\mathbf{F}^{(k+1)}\leftarrow\text{Proj}_{\geq 0}(\mathbf{F}^{(k+1)})\) endif.
6: Obtain \(\mathbf{Z}^{(k+1)}\) by (16).
7:endwhile
8:return Optimal label transport \(\mathbf{F}^{*}\) and model parameters \(\theta^{*}\).
```

**Algorithm 1** Learning a GNN by solving (12)

Denote the objective function in (13) as \(L(\theta,\mathbf{F},\mathbf{Z})\). In the Bregman ADMM framework, we can optimize \(\theta\), \(\mathbf{F}\), and \(\mathbf{Z}\) iteratively through alternating optimization. In the \(k\)-th iteration, we update the three variables via solving the following three subproblems:

\[\begin{split}\theta^{(k+1)}=&\,\text{arg min}_{\theta} \,L(\theta,\mathbf{F}^{(k)},\mathbf{Z}^{(k)})\\ =&\,\text{arg min}_{\theta}\langle\mathbf{Z}^{(k)},\ g _{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta)\rangle+\lambda B_{\phi}(g_{\mathcal{ V}_{L}}(\mathbf{X},\mathbf{A};\theta)+\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}^{(k)},\ \mathbf{Y}_{\mathcal{V}_{L}}).\end{split} \tag{14}\]
\[\mathbf{F}^{(k+1)}= \arg\,\min_{\mathbf{F}\in\mathcal{U}^{|\xi|\times C}}L(\theta^{(k+1)},\bm {F},\mathbf{Z}^{(k)}) \tag{15}\] \[= \arg\,\min_{\mathbf{F}\in\mathcal{U}^{|\xi|\times C}}\langle\mathbf{Z}^{(k )},\ \mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}\rangle+\lambda B_{\phi}(g_{\mathcal{V}_{L}}(\mathbf{X },\mathbf{A};\theta^{(k+1)})+\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F},\ \mathbf{Y}_{\mathcal{V}_{L}}).\]

\[\mathbf{Z}^{(k+1)}=\mathbf{Z}^{(k)}+\lambda(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta^{ (k+1)})+\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}^{(k+1)}-\mathbf{Y}_{\mathcal{V}_{L}}). \tag{16}\]

We can find that (14) is a unconstrained optimization problem, so we can update the model parameter \(\theta\) by gradient descent. Similarly, we can solve (15) and update the flow matrix \(\mathbf{F}\) by gradient descent or projected gradient descent, depending on whether the graph is undirected or not. Finally, the update of the dual variable \(\mathbf{Z}\) can be achieved in a closed form, as shown in (16). The Bregman ADMM algorithm solves the original problem in (12) rather than a relaxed version. In theory, with the increase of iterations, we can obtain the optimal variables that satisfy the equality constraint. Algorithm 2 shows the scheme of the Bregman ADMM-based solver.

#### 3.3.3 Optional Edge Weight Prediction

Some GNNs, e.g., GCN-LPA (Wang and Leskovec, 2021) and GAT (Velickovic et al., 2018), model the adjacency matrix of graph as learnable parameters. Inspired by these models, we can optionally introduce an edge weight predictor and parameterize the adjacency matrix based on the flow matrix, as shown in Figure 1. In particular, given \(\mathbf{F}\), we can apply a multi-layer perceptron (MLP) to embed it to an edge weight vector and then obtain a weighted adjacency matrix. Accordingly, the learning problem becomes

\[\min_{\theta,\xi}\,\min_{\mathbf{F}\in\Omega_{C}(\mathbf{S}_{\mathcal{V}_{L}},g_{ \mathcal{V}_{L}}(\mathbf{X},\mathbf{A}(\mathbf{F};\xi);\theta),\mathbf{Y}_{\mathcal{V}_{L}}) }\|\text{diag}(\mathbf{w})\mathbf{F}\|_{1}, \tag{17}\]

where \(\mathbf{A}(\mathbf{F};\xi)\) represents the adjacency matrix determined by the label transportation \(\mathbf{F}\), and \(\xi\) represents the parameters of the MLP. Taking the learning of \(\xi\) into account, we can modify the above two solvers slightly and make them applicable for solving (17).

### Connections to Traditional Methods

As discussed in Section 3.3.1 and shown in (12), the Bregman divergence \(B_{\phi}\) can be implemented as the commonly-used loss function \(\psi\) in (2) (e.g., the KL divergence or the least-square loss). Table 1 compares the typical setting of traditional learning methods and that of our QW loss-based method in different node-level tasks. Essentially, the traditional learning method in (2) can be viewed as a special case of our QW loss-based learning method. In particular, when setting \(\mathbf{F}=\mathbf{0}_{|\mathcal{F}|\times C}\) and \(B_{\phi}=\psi\), the objective function in (12) degrades to the objective function in (2), which treats each node independently. Similarly, when further setting the dual variable \(\mathbf{Z}=\mathbf{0}_{|\mathcal{V}_{L}|\times C}\), the objective function in (13) degrades to the objective function in (2) as well. In theory, we demonstrate that our QW loss-based learning method can fit training data better than the traditional method does, as shown in the following theorem.

\begin{table}
\begin{tabular}{c|c|c|c} \hline \hline Method & Setting & Node Classification & Node Regression \\ \hline \multirow{2}{*}{\begin{tabular}{c} Apply the \\ loss in (2) \\ \end{tabular} } & \(\psi\) & Cross-entropy or KL & Least-square \\ \cline{2-3}  & Predicted \(\mathbf{y}_{v}\) & \(g_{v}(\mathbf{X},\mathbf{A};\theta)\), \(\forall v\in\mathcal{V}\setminus\mathcal{V}_{L}\) \\ \hline \multirow{2}{*}{
\begin{tabular}{c} Apply the \\ QW loss \\ \end{tabular} } & \(\phi\) & Entropy & \(\frac{1}{2}\|\cdot\|_{2}^{2}\) \\ \cline{2-3}  & \(B_{\phi}(=\psi)\) & KL & Least-square \\ \cline{2-3}  & Predicted \(\mathbf{y}_{v}\) & \(g_{v}(\mathbf{X},\mathbf{A};\theta)+\mathbf{S}_{v}\mathbf{F}\), \(\forall v\in\mathcal{V}\setminus\mathcal{V}_{L}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Comparison between traditional methods and ours
**Theorem 3**.: _Let \(\{\theta^{\star},\mathbf{F}^{\star},\mathbf{Z}^{\star}\}\) be the global optimal solution of (13), \(\{\theta^{\dagger},\mathbf{F}^{\dagger}\}\) be the global optimal solution of (12), and \(\theta^{\dagger}\) be the global optimal solution of (13), \(\{\theta^{\dagger},\mathbf{F}^{\dagger}\}\) be the global optimal solution of (12), and \(\theta^{\dagger}\) be the global optimal solution of \(\min_{\theta}B_{\phi}(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta),\ \mathbf{Y}_{\mathcal{V}_{L}})\), we have_

\[B_{\phi}(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta^{\star})+\mathbf{S}_{\mathcal{V} _{L}}\mathbf{F}^{\star},\ \mathbf{Y}_{\mathcal{V}_{L}})\leq B_{\phi}(g_{\mathcal{V}_{L}}(\mathbf{X},\mathbf{A}; \theta^{\dagger})+\mathbf{S}_{\mathcal{V}_{L}}\mathbf{F}^{\dagger},\ \mathbf{Y}_{\mathcal{V}_{L}})\leq B_{\phi}(g_{ \mathcal{V}_{L}}(\mathbf{X},\mathbf{A};\theta^{\dagger}),\ \mathbf{Y}_{\mathcal{V}_{L}})\]

Proof.: The proof is straightforward -- \(\{\theta^{\ddagger},\mathbf{0}_{|\mathcal{E}|\times C}\}\) is a feasible solution of (12), so the corresponding objective is equal to or larger than that obtained by \(\{\theta^{\dagger},\mathbf{F}^{\dagger}\}\). Similarly, \(\{\theta^{\dagger},\mathbf{F}^{\dagger},\mathbf{0}_{|\mathcal{V}_{L}|\times C}\}\) is a feasible solution of (13), so the corresponding objective is equal to or larger than that obtained by \(\{\theta^{\star},\mathbf{F}^{\star},\mathbf{Z}^{\star}\}\). 

**Remark.** It should be noted that although the Bregman ADMM-based solver can fit training data better in theory, in the cases with distribution shifting or out-of-distribution issues, it has a higher risk of over-fitting. Therefore, in practice, we can select one of the above two solvers to optimize the GW loss, depending on their performance. In the following experimental section, we will further compare these two solvers in details.

### A New Transductive Prediction Paradigm

As shown in Table 1, given the learned model \(\theta^{*}\) and the optimal label transport \(\mathbf{F}^{*}\), we predict node labels in a new transductive prediction paradigm. For \(v\in\mathcal{V}\setminus\mathcal{V}_{L}\), we predict its label as

\[\tilde{\mathbf{y}}_{v}:=g_{v}(\mathbf{X},\mathbf{A};\theta^{*})+\mathbf{S}_{v}\mathbf{F}^{*}, \tag{18}\]

which combines the estimated label from the learned GNN and the residual component from the optimal label transport.

It should be noted that some attempts have been made to incorporate label propagation algorithms (LPAs) (Zhu and Goldberg, 2022) into GNNs, e.g., the GCN-LPA in (Wang and Leskovec, 2021) and the FDiff-Scale in (Huang et al., 2021). The PTA in (Dong et al., 2021) demonstrates that learning a decoupled GNN is equivalent to implementing a label propagation algorithm. These methods leverage LPAs to regularize the learning of GNNs. However, in the prediction phase, they abandon the training labels and rely only on the GNNs to predict node labels, as shown in Table 1. Unlike these methods, our QW loss-based method achieves a new kind of label propagation with the help of computational optimal transport, saving the training label information in the optimal label transport and applying it explicitly in the prediction phase.

## 4 Experiments

To demonstrate the effectiveness of our QW loss-based learning method, we apply it to learn GNNs with various architectures and test the learned GNNs in different node-level prediction tasks. We compare our learning method with the traditional one in (2) on their model performance and computational efficiency. For our method, we conduct a series of analytic experiments to show its robustness to hyperparameter settings and label insufficiency. All the experiments are conducted on a machine with three NVIDIA A40 GPUs, and the code is implemented based on PyTorch.

### Implementation Details

#### 4.1.1 Datasets

The datasets we considered consist of five homophilic graphs (i.e., **Cora**, **Citeseer**, **Pubmed**(Sen et al., 2008; Yang et al., 2016), **Computers**, and **Photo**(McAuley et al., 2015)) and five heterophilic
graphs (i.e., **Chameleon**, **Squirrel**(Rozemberczki et al., 2021), **Actor**, **Texas**, and **Cornell**(Pei et al., 2020)), respectively. Following the work in (Wang and Leskovec, 2021), we categorize the graphs according to the percentage of the edges connecting the nodes of the same class (i.e., the intra-edge rate). The basic information of these graphs is shown in Table 2. Additionally, a large **arXiv-year** graph (Lim et al., 2021) is applied to demonstrate the efficiency of our method. The adjacency matrix of each graph is binary, so the edge weights \(\mathbf{w}=\mathbf{1}_{|\mathcal{E}|}\).

#### 4.1.2 GNN Architectures

In the following experiments, the models we considered include \(i)\) the representative spatial GNNs, i.e., **GCN**(Kipf and Welling, 2017), **GAT**(Velickovic et al., 2018), **GIN**(Xu et al., 2019), **GraphSAGE**(Hamilton et al., 2017), and **GCN-LPA**(Wang and Leskovec, 2021) that combines the GCN with the label propagation algorithm; and \(ii)\) state-of-the-art spectral GNNs, i.e., **APPNP**(Gasteiger et al., 2019), **BernNet**(He et al., 2021), and **ChebNetII**(He et al., 2022). For a fair comparison, we set the architectures of the GNNs based on the code provided by (He et al., 2022) and configure the algorithmic hyperparameters by grid search. More details of the hyperparameter settings are in Appendix.

#### 4.1.3 Learning Tasks and Evaluation Measurements

For each graph, their nodes belong to different classes. Therefore, we first learn different GNNs to solve the node-level classification tasks defined on the above graphs. By default, the split ratio of each graph's nodes is 60% for training, 20% for validation, and 20% for testing, respectively. The GNNs are learned by \(i)\) the traditional learning method in (2)1 and \(ii)\) minimizing the proposed QW loss in (11), respectively. When implementing the QW loss, we apply either Algorithm 1 or 2, depending on their performance. Additionally, to demonstrate the usefulness of our QW loss in node-level regression tasks, we treat the node labels as one-hot vectors and fit them by minimizing the mean squared error (MSE), in which the \(\psi\) in (2) and the corresponding Bregman divergence

\begin{table}
\begin{tabular}{c|c c c c c|c c c c c|c} \hline \hline \multirow{2}{*}{Model Method} & \multicolumn{6}{c|}{Homophilic graphs} & \multicolumn{6}{c|}{Heterophilic graphs} \\ \cline{2-13}  & Cora & Citeseer & Pubmed & Computers & Photo & Squirrel & Chameleon & Actor & Texas & Cornell \\ \hline \#Nodes (\(|\mathcal{V}|\)) & 2,708 & 3,327 & 19,717 & 13,752 & 5,201 & 7,650 & 2,277 & 7,600 & 183 & 183 & Overall \\ \#Features (\(D\)) & 1,433 & 3,703 & 500 & 767 & 754 & 2,089 & 2,325 & 932 & 1,703 & 1,703 & Improve \\ \#Edaps (\(|\mathcal{E}|\)) & 5,278 & 4,552 & 44,324 & 25,861 & 19,081 & 19,358 & 31,371 & 26,659 & 279 & 277 & Improve \\ Intra-edge rate & 81.0\% & 73.6\% & 80.2\% & 77.7\% & 82.7\% & 22.2\% & 23.0\% & 21.8\% & 6.1\% & 12.3\% & \\ \#Classes (\(C\)) & 7 & 6 & 5 & 10 & 8 & 5 & 5 & 5 & 5 & 5 & 5 \\ \hline \multirow{2}{*}{GCN} & (2) & 87.44\({}_{10.39}\) & 79.93\({}_{10.84}\) & 86.93\({}_{10.29}\) & 88.42\({}_{10.43}\) & 93.24\({}_{10.43}\) & 65.55\({}_{11.51}\) & 63.57\({}_{11.06}\) & 34.00\({}_{12.8}\) & 77.21\({}_{11.28}\) & 61.91\({}_{15.11}\) & — \\  & (2)+LPM & 86.34\({}_{10.4}\) & 78.51\({}_{11.22}\) & 84.72\({}_{10.70}\) & 82.48\({}_{10.80}\) & 88.10\({}_{11.31}\) & 44.81\({}_{11.81}\) & 60.90\({}_{10.53}\) & 32.43\({}_{11.50}\) & 78.09\({}_{14.74}\) & 67.27\({}_{23.50}\) & -1.36 \\  & QW & **87.88\({}_{10.78}\)** & **81.36\({}_{10.41}\)** & **87.89\({}_{10.40}\)** & **89.20\({}_{10.41}\)** & **93.81\({}_{10.56}\)** & **25.62\({}_{10.49}\)** & **68.10\({}_{11.01}\)** & **38.09\({}_{10.59}\)** & **84.10\({}_{10.29}\)** & **84.26\({}_{10.38}\)** & **4.26\({}_{10.38}\)** & **4.81\({}_{10.29}\)** & **84.26\({}_{10.38}\)** \\ \hline \multirow{2}{*}{ GAT} & (2) & **89.20\({}_{10.78}\)** & **87.07\({}_{10.78}\)** & **87.42\({}_{10.33}\)** & 90.80\({}_{10.36}\)** & **89.35\({}_{10.29}\)** & **84.35\({}_{10.29}\)** & **85.22\({}_{10.41}\)** & 64.31\({}_{11.21}\) & **35.68\({}_{10.80}\)** & 80.00\({}_{11.69}\)** & **89.21\({}_{13.13}\)** & **11.41\({}_{11.23}\)** & **11.41** \\  & QW & 89.11\({}_{11.01}\) & 89.10\({}_{10.84}\) & **88.50\({}_{10.29}\)** & **90.01\({}_{10.28}\)** & **94.65\({}_{10.28}\)** & **95.64\({}_{10.28}\)** & **95.53\({}_{11.39}\)** & **67.35\({}_{11.42}\)** & **35.66\({}_{11.32}\)** & **33.66\({}_{11.23}\)** & **70.21\({}_{11.23}\)** & **11.14** \\ \hline \multirow{2}{*}{GIN} & (2) & 86.22\({}_{10.20}\) & **76.18\({}_{10.78}\)** & **87.87\({}_{10.23}\)** & 80.87\({}_{11.43}\) & 89.83\({}_{10.72}\) & 39.11\({}_{11.23}\) & 64.29\({}_{11.51}\) & **32.37\({}_{11.56}\)** & **27.79\({}_{11.42}\)** & 62.55\({}_{10.54}\)** & **-** \\  & QW & **86.24\({}_{10.00}\)** & **76.13\({}_{11.00}\)** & **87.53\({}_{10.34}\)** & **89.28\({}_{10.29}\)** & **92.84\({}_{10.26}\)** & **92.60\({}_{10.46}\)** & **65.98\({}_{10.26}\)** & **73.26\({}_{11.32}\)** & **32.32\({}_{11.32}\)** & **77.54\({}_{10.42}\)** & **64.40\({}_{10.42}\)** & **-** \\ \hline \multirow{2}{*}{GraphSAGE} & (2) & 88.24\({}_{10.49}\) & 79.51\({}_{10.80}\) & 88.53\({}_{10.34}\) & 89.87\({}_{11.01}\) & 95.50\({}_{10.42}\) & 43.79\({}_{10.90}\) & 63.26\({}_{10.99}\) & **89.90\({}_{10.90}\)** & **80.30\({}_{10.34}\)** & **86.32\({}_{10.38}\)** & **-** \\  & QW & 87.59\({}_{10.77}\) & **80.52\({}_{10.68}\)** & **88.61\({}_{10.39}\)** & **90.17
\(B_{\phi}\) are set to be the least-square loss. For each method, we perform 10 runs with different seeds and record the learning results' mean and standard deviation.

### Numerical Comparison and Visualization

#### 4.2.1 Node Classification and Regression

Table 2 shows the node classification results on the ten graphs,2 whose last column records the overall improvements caused by our QW loss compared to other learning methods. The experimental results demonstrate the usefulness of our QW loss-based learning method -- for each model, applying our QW loss helps to improve learning results in most situations and leads to consistent overall improvements. In particular, for the state-of-the-art spectral GNNs like BernNet (He et al., 2021) and ChebNetII (He et al., 2022), learning with our QW loss can improve their overall performance on both homophilic and heterophilic graphs consistently, resulting in the best performance in this experiment. For the simple GCN model (Kipf and Welling, 2017), learning with our QW loss improves its performance significantly and reduces the gap between its classification accuracy and that of the state-of-the-art models (Gasteiger et al., 2019; He et al., 2021, 2022), especially heterophilic graphs. Note that, when learning GCN, our QW loss works better than the traditional method regularized by the label propagation (i.e., GCN-LPA (Wang and Leskovec, 2021)) because our learning method can leverage the training label information in both learning and prediction phases. GCN-LPA improves GCN when learning on heterophilic graphs, but surprisingly, leads to performance degradation on homophilic graphs. Additionally, we also fit the one-hot labels by minimizing the MSE. As shown in Table 3, minimizing the QW loss leads to lower MSE results, which demonstrates the usefulness of the QW loss in node-level regression tasks.

Footnote 2: In Table 2, we bold the best learning result for each graph. Learning GCN by “(2)+LPA” means implementing GCN-LPA (Wang and Leskovec, 2021).

#### 4.2.2 Computational Efficiency and Scalability

In theory, the computational complexity of our QW loss is linear with the number of edges. When implementing the loss as (13) and solving ti by Algorithm 2, its complexity is also linear with the number of inner iterations \(J\). Figure 2 shows the runtime comparisons for the QW loss-based learning methods and the traditional method on two datasets. We can find that minimizing the QW loss by Algorithm 1 or Algorithm 2 with \(J=1\) merely increases the training time slightly compared to the traditional method. Empirically, setting \(J\leq 5\) can leads to promising learning results, as shown in Tables 2 and 3. In other words, the computational cost of applying the QW loss is tolerable considering the significant performance improvements it achieved. Additionally, we apply our QW loss to large-scale graphs and test its scalability. As shown in Table 4, we implement the QW loss based on Algorithm 1 (i.e., solving (12)) and apply it to the node classification task

\begin{table}
\begin{tabular}{c c|c c|c c} \hline \hline \multirow{2}{*}{Model} & \multirow{2}{*}{Method} & \multicolumn{2}{c|}{Homophilic graphs} & \multicolumn{2}{c}{Heterophilic graphs} \\  & & Computers & Photo & Actor & Cornell \\ \hline \multirow{2}{*}{GIN} & (2) & \(0.0605_{\pm 0.0018}\) & \(0.0459_{\pm 0.0044}\) & \(0.1570_{\pm 0.0014}\) & \(0.1609_{\pm 0.0359}\) \\  & QW & \(\mathbf{0.0244_{\pm 0.0028}}\) & \(\mathbf{0.0203_{\pm 0.0012}}\) & \(\mathbf{0.1564_{\pm 0.0012}}\) & \(\mathbf{0.1524_{\pm 0.0043}}\) \\ \hline \multirow{2}{*}{BernNet} & (2) & \(0.0871_{\pm 0.0002}\) & \(0.0488_{\pm 0.0009}\) & \(\mathbf{0.1661_{\pm 0.0020}}\) & \(0.0989_{\pm 0.0076}\) \\  & QW & \(\mathbf{0.0364_{\pm 0.0038}}\) & \(\mathbf{0.0297_{\pm 0.0014}}\) & \(0.1671_{\pm 0.0008}\) & \(\mathbf{0.0753_{\pm 0.0024}}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Comparisons on node regression error (MSE).
in the large-scale arXiv-year graph. The result shows that our QW loss is applicable to the graphs with millions of edges on a single GPU and improves the model performance.

#### 4.2.3 Distribution of Optimal Label Transport

Figure 3 visualizes the histograms of the optimal label transport \(\mathbf{F}^{*}\) learned for two representative GNNs (i.e., GCN (Kipf and Welling 2017) and ChebNetII (He et al. 2022)) on four graphs (i.e., the homophilic graphs "Computers" and "Photo" and the heterophilic graphs "Squirrel" and "Chameleon"). We can find that when learning on homophilic graphs, the elements of the optimal label transport obey the zero-mean Laplacian distribution. It is reasonable from the perspective of optimization -- the term \(\|\text{diag}(\mathbf{w})\mathbf{F}\|_{1}\) can be explained as a Laplacian prior imposed on \(\mathbf{F}\)'s element. Additionally, we can find that the distribution corresponding to GCN has larger variance than that corresponding to ChebNetII, which implies that the \(\mathbf{F}^{*}\) of GCN has more non-zero elements and thus has more significant impacts on label prediction. The numerical results in Table 2 can also verify this claim -- in most situations, the performance improvements caused by the optimal label transport is significant for GCN but slight for ChebNetII. For heterophilic graphs, learning GCN still leads to Laplacian distributed label transport. However, the distributions corresponding to ChebNetII are diverse -- the distribution for Squirrel is long-tailed while that for Chameleon is still Laplacian.
### Analytic Experiments

#### 4.3.1 Robustness to Label Insufficiency Issue

Our QW loss-based learning method considers the label transport on graphs, whose feasible domain is determined by the observed training labels. The more labels we observed, the smaller the feasible domain is. To demonstrate the robustness of our method to the label insufficiency issue, we evaluate the performance of our method given different amounts of training labels. We train ChebNetII (He et al., 2022) on four graphs by traditional method and our method, respectively. For each graph, we use \(K\%\) nodes' labels to train the ChebNetII, where \(K\in[1,60]\), and apply \(20\%\) nodes for validation and \(20\%\) nodes for testing, respectively, as the above default setting does. Figure 4 shows that our QW loss-based learning method can achieve encouraging performance even when only \(20\%\) nodes or fewer are labeled. Additionally, the methods are robust to the selection of solver -- we can minimize the QW loss based on (12) or (13), leading to comparable results and outperforming the traditional loss consistently.

#### 4.3.2 Impacts of Adjusting Edge Weights

As shown in (17), we can train an MLP to predict edge weights based on the optimal label transport. The ablation study in Table 5 quantitatively show the impacts of adjusting edge weights on the learning results. We can find that for ChebNetII, the two settings provide us with comparable learning results. For GCN, however, learning the model with adjusted edge weights suffers from performance degradation on homophilic graphs while leads to significant improvements on heterophilic graphs. Empirically, it seems that adjusting edge weights based on label transportation helps to improve the learning of simple GNN models on heterophilic graphs.

\begin{table}
\begin{tabular}{c c|c c|c c} \hline \hline \multicolumn{2}{c|}{Model} & \multicolumn{2}{c|}{Homophilic} & \multicolumn{2}{c}{Heterophilic} \\  & & Computers & Photo & Actor & Cornell \\ \hline \multirow{2}{*}{GCN} & \(\mathbf{A}\) & \(\mathbf{88.39_{\pm 0.55}}\) & \(\mathbf{93.80_{\pm 0.37}}\) & \(30.14_{\pm 0.80}\) & \(60.64_{\pm 4.26}\) \\  & \(\mathbf{A}(\mathbf{F};\xi)\) & \(84.35_{\pm 0.46}\) & \(91.79_{\pm 0.21}\) & \(\mathbf{38.09_{\pm 0.50}}\) & \(\mathbf{84.26_{\pm 2.98}}\) \\ \hline \multirow{2}{*}{ChebNetII} & \(\mathbf{A}\) & \(\mathbf{89.52_{\pm 0.54}}\) & \(\mathbf{94.84_{\pm 0.37}}\) & \(\mathbf{41.37_{\pm 0.67}}\) & \(86.38_{\pm 3.19}\) \\  & \(\mathbf{A}(\mathbf{F};\xi)\) & \(89.41_{\pm 0.41}\) & \(94.79_{\pm 0.45}\) & \(40.74_{\pm 0.80}\) & \(\mathbf{86.60_{\pm 2.98}}\) \\ \hline \hline \end{tabular}
\end{table}
Table 5: Impacts of adjusting edge weights on node classification accuracy (%) when applying the QW loss.

Figure 4: Illustrations of the learning methods’ performance given different amounts of labeled nodes.
#### 4.3.3 Robustness to Hyperparameters

The weight of the Bregman divergence term, i.e., \(\lambda\), is the key hyperparameter impacting the performance of our learning method. Empirically, when \(\lambda\) is too small, the Bregman divergence between the observed labels and their predictions becomes ignorable. Accordingly, the regularizer may be too weak to supervise GNNs' learning properly. On the contrary, when \(\lambda\) is too large, the regularizer becomes dominant in the learning objective, and the impact of the label transport becomes weak in both the learning and prediction phases. As a result, it may perform similarly to the traditional method when using a large \(\lambda\). We test the robustness of our method to \(\lambda\) and show representative results in Figure 5. In particular, our QW loss-based method trains ChebNetII (He et al., 2022) on four graphs. Both Algorithm 1 for (12) and Algorithm 2 for (13) are tested. The \(\lambda\) is set in the range from \(10^{-2}\) to \(10^{3}\). For homophilic graphs, our learning method achieves stable performance when \(\lambda\geq 10\). When \(\lambda<10\), the learning results degrade significantly because of inadequate supervision. Our learning method often obtains the best learning result for heterophilic graphs when \(\lambda\in[1,10]\). These experimental results show that our method is robust to the setting of \(\lambda\), and we can set \(\lambda\) in a wide range to obtain relatively stable performance.

## 5 Conclusion

We have proposed the Quasi-Wasserstein loss for learning graph neural networks. This loss matches well with the non-i.i.d. property of graph-structured data, providing a new strategy to leverage observed node labels in both training and testing phases. Applying the QW loss to learn GNNs improves their performance in various node-level prediction tasks. In the future, we would like to explore the impacts of the optimal label transport on the generalization power of GNNs in theory. Moreover, we plan to modify the QW loss further, developing a new optimization strategy to accelerate its computation.
