Improving Graph Neural Networks by Learning Continuous Edge Directions

## Abstract

Abstract Graph Neural Networks (GNNs) traditionally employ a message-passing mechanism that resembles diffusion over undirected graphs, which often leads to homogenization of node features and reduced discriminative power in tasks such as node classification. Our key insight for addressing this limitation is to assign fuzzy edge directions—that can vary continuously from node i 𝑖 i italic_i pointing to node j 𝑗 j italic_j to vice versa—to the edges of a graph so that features can preferentially flow in one direction between nodes to enable long-range information transmission across the graph. We also introduce a novel complex-valued Laplacian for directed graphs with fuzzy edges where the real and imaginary parts represent information flow in opposite directions. Using this Laplacian, we propose a general framework, called Continuous Edge Direction (CoED) GNN, for learning on graphs with fuzzy edges and prove its expressivity limits using a generalization of the Weisfeiler-Leman (WL) graph isomorphism test for directed graphs with fuzzy edges. Our architecture aggregates neighbor features scaled by the learned edge directions and processes the aggregated messages from in-neighbors and out-neighbors separately alongside the self-features of the nodes. Since continuous edge directions are differentiable, they can be learned jointly with the GNN weights via gradient-based optimization. CoED GNN is particularly well-suited for graph ensemble data where the graph structure remains fixed but multiple realizations of node features are available, such as in gene regulatory networks, web connectivity graphs, and power grids. We demonstrate through extensive experiments on both synthetic and real datasets that learning continuous edge directions significantly improves performance both for undirected and directed graphs compared with existing methods.

## 1 Introduction

Graph Neural Networks (GNNs) have emerged as a powerful tool for learning from data that is structured as graphs, with applications ranging from social network analysis to molecular chemistry (Kipf & Welling, ; Zhou et al., ; Gilmer et al., ). GNNs typically employ a message passing mechanism where nodes aggregate and then transform feature information from their neighbors at each layer, enabling them to learn node representations that capture both local and global graph structures. When the graph is undirected, the aggregation of node features mimics a diffusion process. Each node’s representation becomes the averaged features of its immediate neighbors, leading to a homogenization of information across the graph. As depth increases, this diffusion of information culminates in a uniform state where node representations converge towards a constant value across all the nodes, which severely limits the discriminative power of GNNs, especially in tasks such as node classification (Rusch et al., ; Oono & Suzuki, ; Cai & Wang, ; Li et al., ; Keriven, ; Chen et al., ; Wu et al., ; ).

Figure: Figure 1: (a) When edges are undirected information diffuses across the graph and long range transmission of information between nodes 1 and 2 is not possible. (b) Once the optimal edge directions are learned, information can flow directly from node 1 to node 2.
Refer to caption: x1.png

Our key insight for improving the performance of GNNs is to alter the nature of information transmission between nodes from diffusion to flow. To do so, we add directions to the edges of a graph so that features can be propagated from node $v_{i}$ to its neighbor node $v_{j}$ without a reciprocal propagation of information from node $v_{j}$ to node $v_{i}$. Unlike diffusion, where information uniformly spreads across available paths, flow is directional and preserves the propagation of information across longer distances within a graph, as illustrated in Figure 1. In general, the optimal information propagation could require edges whose directions fall anywhere in the continuum of node $v_{i}$ pointing to node $v_{j}$ to vice versa.

To capture such continuous edge directions, we propose a concept of ‘fuzzy edges,’ where the direction of an edge between any two nodes $v_{i}$ and $v_{j}$ is not a discrete but a continuous value. An edge’s orientation can range continuously—from exclusively pointing from node $v_{i}$ to node $v_{j}$, through a fully bidirectional state, to exclusively pointing from node $v_{j}$ to node $v_{i}$. Therefore, ‘fuzzy’ direction essentially controls the relative amount of information flow from node $v_{i}$ to node $v_{j}$ and the reciprocal flow from node $v_{j}$ to node $v_{i}$. To effectively model this directional flexibility, we introduce a complex-valued graph Laplacian called a fuzzy Laplacian. In this framework, the real part of the $ij$-th entry in the fuzzy Laplacian matrix quantifies the degree of information transmission from node $v_{j}$ to node $v_{i}$, while the imaginary part measures the flow from node $v_{i}$ to node $v_{j}$.

Next, we introduce the Continuous Edge Direction (CoED) GNN architecture. At each layer, a node’s neighbors’ features are scaled by the directions of their connecting edges and aggregated. This aggregation is performed separately for incoming and outgoing edges, following Rossi et al. (), resulting in distinct features for incoming and outgoing messages. Practically, this is implemented by applying the fuzzy Laplacian to the node features, where the real and imaginary parts correspond to the features aggregated from incoming neighbors and outgoing neighbors, respectively. These aggregated features are then affine transformed using learnable weights and combined with the node’s own transformed features. A nonlinear activation function is applied to obtain the updated node features. This process is repeated for each layer. The continuous edge directions have the added benefit that they are differentiable. During training, both the edge directions and weight matrices are learned simultaneously using gradient-based optimization to improve the learning objective.

Importantly, our approach is fundamentally different from methods such as Graph Attentions Network (GAT) (Veličković et al., ) or graph transformers (Dwivedi & Bresson, ; Rampášek et al., ) that learn attention coefficients to assign weights to each edge of the graph based on the features (and potentially the positional encoding) of the nodes connected by that edge. While such an attention mechanism can capture asymmetric relationships by computing direction-specific attention weights based on node features, they do not learn edge directions as independent parameters. In these models, the attention coefficient from node $v_{i}$ to node $v_{j}$ are functions of the features of
$v_{i}$ and $v_{j}$, and will change if node features change. In contrast, our approach introduces continuous edge directions as learnable parameters that are optimized end-to-end, independent of the node features. As we demonstrate empirically below, by directly learning edge directions, our method goes beyond what can be achieved through an attention mechanism, enabling long-range information flow on graphs and improved performance.

End-to-end learning of edge directions is most effective on graph ensemble data, where the graph structure remains fixed but multiple realizations of node features and targets (such as node labels) are available. This effectiveness arises from the ability to optimize information flow across all edges simultaneously without a need to mask parts of the graph for training and testing. Instead, training and testing splits are based on different feature realizations rather than on subsets of the graph. Graph ensemble data are increasingly common across various domains. For example, in biology, gene-regulatory networks are constant directed graphs where nodes represent genes and edges represent gene-gene interactions, while node features like gene expression levels vary across different cells. Similarly, in web connectivity, the network of websites remains relatively static, but traffic patterns change over time, providing different node feature sets on the same underlying graph. In power grids, the network of electrical components is fixed, while the steady-state operating points of these components vary under different conditions, yielding multiple observations on the same graph. In all these cases, a fixed graph is paired with numerous feature variations. By applying CoED GNN to these scenarios, we demonstrate that learning edge directions significantly improves performance for both directed and undirected graphs. The main contributions of this paper are the following:

- •
We introduce a principled complex-valued graph Laplacian for graphs where edge directions can vary continuously and prove that it is more expressive than existing forms of Laplacians for directed graphs, such as the magnetic Laplacian.
- •
We propose an architecture called Continuous Edge Direction (CoED) GNN, which is a general framework for learning on directed graphs with fuzzy edges. We prove that CoED GNN is as expressive as a weak form of the Weisfeiler-Leman (WL) graph isomorphism test for directed graphs with fuzzy edges.
- •
Using extensive experiments, we empirically show that learning edge directions significantly improves performance by applying CoED GNN to both synthetic and real graph ensemble data.

## 2 Preliminaries

A graph is defined as a pair $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{v_{1},v_{2},\ldots,v_{N}\}$ is a set of $N$ nodes, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is a set of edges connecting pairs of nodes. Each node $v_{i}$ is associated with a feature vector $\mathbf{f}_{i}\in\mathbb{R}^{D}$, where $D$ is the dimensionality of the feature space, and collectively these feature vectors form the node feature matrix $\mathbf{F}\in\mathbb{R}^{N\times D}$. Additionally, each node is assigned a prediction target, such as a class label for classification tasks or a continuous value for regression tasks.

The connectivity of a graph is encoded in an adjacency matrix $\mathbf{A}\in\{0,1\}^{N\times N}$. If there is an undirected edge between $v_{i}$ and $v_{j}$, then both $A_{ij}=1$ and $A_{ji}=1$. For a directed edge, one of $A_{ij}$ or $A_{ji}$ is 1 while the other is 0, specifying the direction of information flow. A directed edge, $A_{ij}=1$ and $A_{ji}=0$, indicates that $v_{j}$ sends information to $v_{i}$. Hence, we refer to $v_{j}$ as the in-neighbor of $v_{i}$ and conversely to $v_{i}$ as the out-neighbor of $v_{j}$. If both $A_{ij}=0$ and $A_{ji}=0$, there is no edge between $v_{i}$ and $v_{j}$. Accordingly, in a directed graph, we define two distinct degree matrices: the in-degree matrix $\mathbf{D_{in}}=\mathrm{diag}(\mathbf{A}\mathbf{1})$ and the out-degree matrix $\mathbf{D_{out}}=\mathrm{diag}(\mathbf{A}^{\top}\mathbf{1})$.

In GNNs, node features $\mathbf{F}$ are processed iteratively through a message-passing mechanism that leverages the structural information of the graph $\mathcal{G}$. This process involves two main steps at each layer $l$:

1. Message Aggregation: For each node $v_{i}$, an aggregated message $\mathbf{m}_{i,\mathcal{N}(i)}^{(l)}$ is computed from the features of its neighbors:

$$ $\mathbf{m}_{i,\mathcal{N}(i)}^{(l)}=\text{AGGREGATE}\left(\{\!\!\{(\mathbf{f}_ {i}^{(l-1)},\mathbf{f}_{j}^{(l-1)})\mid j\in\mathcal{N}(i)\}\!\!\}\right)$ $$

Here, $\mathcal{N}(i)$ denotes the set of nodes $v_{j}$ that are connected to node $v_{i}$ by an edge.

2. Feature Update: The feature vector of node $v_{i}$ is then updated using the aggregated message:

$$ $\mathbf{f}_{i}^{(l)}=\text{UPDATE}\left(\mathbf{f}_{i}^{(l-1)},\mathbf{m}_{i, \mathcal{N}(i)}^{(l)}\right)$ $$

AGGREGATE and UPDATE are functions with learnable parameters, and their specific implementations define different GNN architectures (Gilmer et al., ).

## 3 Formulation of GNN on Directed Graphs with Fuzzy Edges

### 3.1 Continuous edge directions as phase angles

To describe a continuously varying edge direction between node $v_{i}$ and node $v_{j}$, we assign an angle $\theta_{ij}\in[0,\pi/2]$ to the edge connecting $v_{i}$ to $v_{j}$. During aggregation of features from neighbors, features propagated from $v_{j}$ to $v_{i}$ are scaled by a factor of $\cos{\theta_{ij}}$. Conversely, the features that $v_{j}$ receives from $v_{i}$ are scaled by $\sin{\theta_{ij}}$. For example, when $\theta_{ij}=0$, we have a directed edge where $v_{j}$ sends messages to $v_{i}$ but does not receive any messages from $v_{i}$. When $\theta_{ij}=\pi/4$, the edge is undirected and the same scaling is applied to the messages sent and received by $v_{i}$ to and from $v_{j}$, i.e., $\cos{\pi/4}=\sin{\pi/4}=1/\sqrt{2}$. To ensure consistency, we require that the message received by $v_{i}$ from $v_{j}$ should be equivalent to the message sent by $v_{j}$ to $v_{i}$. It follows that $\theta_{ji}=\pi/2-\theta_{ij}$. We define the phase matrix $\Theta\in[0,\pi/2]^{N\times N}$ to describe the directions of all the edges in a graph. $(\Theta)_{ij}$ is only defined if there is an edge connecting nodes $v_{i}$ and $v_{j}$.

### 3.2 Fuzzy Graph Laplacian

To keep our message-passing GNN as expressive as possible, we define a Laplacian matrix that, during the aggregation step, propagates information along directed edges but keeps the aggregated features from in-neighbors and out-neighbors for each node distinct by assigning them to the real and imaginary parts of a complex number, respectively. For a given $\Theta$, we construct the corresponding fuzzy graph Laplacian $\mathbf{L}_{F}$ as follows. The diagonal entries $(\mathbf{L}_{F})_{ii}$ are zero as we cannot define edge directions for self-loops, and off-diagonal entries are either zero or a phase value,

$$ $(\mathbf{L}_{F})_{ij}=\begin{cases}0&\text{if }A_{ij}=A_{ji}=0\\ \exp(i\theta_{ij})&\text{otherwise}\end{cases}$ (1) $$

Since $\theta_{ij}$ and $\theta_{ji}$ are related by $\theta_{ji}=\pi/2-\theta_{ij}$, it follows that $\mathbf{L}_{F}=i\mathbf{L}_{F}^{\dagger}$, where $\dagger$ is the conjugate transpose. $\text{Re}[{\mathbf{L}_{F}}]$ thus encodes all $i\leftarrow j$ edges scaled by $\cos{\theta_{ij}}$ and $\text{Im}[{\mathbf{L}_{F}}]$ all $i\rightarrow j$ edges scaled by $\sin{\theta_{ij}}$. In Appendix D, we show the fuzzy graph Laplacian admits orthogonal eigenvectors with eigenvalues of the form $a+ia$ with $a\in\mathbb{R}$. Therefore, the eigenvectors of our Laplacian provide positional encodings that are informed by the directions of the edges in addition to their connectivities, as shown in Appendix C. In addition, we prove in Appendix E that a message-passing GNN whose aggregation step is done using the fuzzy Laplacian is as expressive as an extension of the Weisfeiler-Leman (WL) graph isomorphism test to directed graphs with fuzzy edges.

An alternative form of Laplacian for directed graphs is the magnetic Laplacian (Zhang et al., ; He et al., ), which is also complex-valued. For directed graphs with fuzzy edges, the magnetic Laplacian is not as expressive as the Laplacian proposed above. We provide a proof in Appendix F. Briefly, the magnetic Laplacian produces linear combinations of in- and out- neighbor messages at each node as both the real and imaginary parts of the aggregated features. In principle, GNNs should be able to disentangle these linear combinations to recover the in- and out- neighbor messages. However, the linear combinations depend on the local neighborhood of each node which is distinct from one node to another, whereas GNN parameters are shared across all nodes. Therefore, in general, a GNN using the magnetic Laplacian loses the ability to disentangle the in- and out- neighbor messages at each node and thus has lower expressivity. The Laplacian proposed above does not suffer from this limitation since by construction the real and imaginary values directly correspond to the in- and out- neighbor aggreagated messages, respectively.

### 3.3 Model Architecture: Continuous Edge Direction (CoED) GNN

To ensure maximum expressivity, a message-passing mechanism on a directed graph should, for each node, separately aggregate the features of the in-neighbors and the out-neighbors, and independently process each of the aggregated features and the self-features to obtain an updated feature for each node. To this end, we define in- and out- edge weight matrices as $\mathbf{A_{\leftarrow}}=\text{Re}[{\mathbf{L}_{F}}]$ and $\mathbf{A_{\rightarrow}}=\text{Im}[{\mathbf{L}_{F}}]$, respectively. We compute the in- and out- degree matrices as $\mathbf{D_{\leftarrow}}=\mathrm{diag}\left(\mathbf{A_{\leftarrow}}\mathbf{1}\right)$ and $\mathbf{D_{\rightarrow}}=\mathrm{diag}\left(\mathbf{A_{\rightarrow}}\mathbf{1}\right)$. Following Rossi et al. (), but extending it to graphs with continuous edge directions, we define in- and out- fuzzy propagation matrices as

$$ $\begin{split}\mathbf{P}_{\leftarrow}=\mathbf{D}^{-1/2}_{\leftarrow}\mathbf{A}_ {\leftarrow}\mathbf{D}^{-1/2}_{\rightarrow}\\ \mathbf{P}_{\rightarrow}=\mathbf{D}^{-1/2}_{\rightarrow}\mathbf{A}_{ \rightarrow}\mathbf{D}^{-1/2}_{\leftarrow}\end{split}$ (2) $$

Using these matrices, we compute in- and out- messages at layer $l$ as

$$ $\begin{split}\mathbf{m}^{(l)}_{\leftarrow}=\mathbf{P}_{\leftarrow}\mathbf{F}^{ (l-1)}\\ \mathbf{m}^{(l)}_{\rightarrow}=\mathbf{P}_{\rightarrow}\mathbf{F}^{(l-1)}\end{split}$ (3) $$

which defines the AGGREGATE function.

Since self-loops are omitted from $\mathbf{L}_{F}$, we include current node features along with the two directional messages in UPDATE function and update node features as,

$$ $\mathbf{F}^{(l)}=\sigma\big{(}\mathbf{F}^{(l-1)}\mathbf{W}_{\text{self}}^{(l)} +\mathbf{m}^{(l)}_{\leftarrow}\mathbf{W}_{\leftarrow}^{(l)}+\mathbf{m}^{(l)}_{ \rightarrow}\mathbf{W}_{\rightarrow}^{(l)}+\mathbf{B}^{(l)}\big{)}$ (4) $$

where $\sigma$ is an activation function, and $\mathbf{W}_{\text{self}/\leftarrow/\rightarrow}^{(l)}$ and $\mathbf{B}^{(l)}$ are self/in/out weight matrices and a bias matrix, respectively.

The features at the final layer are then transformed using a linear layer to obtain the output for a specific learning task, e.g. node classification or node regression. We use end-to-end gradient-based optimization to iteratively update both the phase matrix $\Theta$ and the GNN parameters $\mathbf{W}_{\text{self}/\leftarrow/\rightarrow}^{(l)}$ and $\mathbf{B}^{(l)}$ at each layer, as illustrated in Figure 2. We allow for the option to learn a different set of edge directions at each layer, $\Theta^{(l)}$, just as we have distinct GNN parameters at each layer.

Figure: Figure 2: Schematic of training with a graph ensemble data. The input graph is undirected (left box). The graph ensemble data contains multiple realizations of node features and corresponding target values, either at the node, edge, or graph level. The phase angle formulation allows continuous edge directions to be optimized alongside the GNN parameters in an end-to-end manner (middle box). The learned edge directions (right box) enable long range information transmission across the graph.
Refer to caption: x2.png

## 4 Related Work

The issue of feature homogenization in GNNs, known as the oversmoothing problem, has been a significant concern. Early studies identified the low-pass filtering effect of GNNs (Defferrard et al., ; Wu et al., ), linking it to oversmoothing and loss of discriminative power (Li et al., ; Oono & Suzuki, ). Proposed solutions include regularization techniques like edge dropout (Rong et al., ), feature masking (Hasanzadeh et al., ), layer normalization (Zhao & Akoglu, ), incorporating signed edges (Derr et al., ), adding residual connections (Chen et al., ), gradient gating (Rusch et al., ), and constraining the Dirichlet energy (Zhou et al., ). Dynamical systems approaches have also been explored, modifying message passing via nonlinear heat equations (Eliasof et al., ), coupled oscillators (Rusch et al., ), and interacting particle systems (Wang et al., ; Di Giovanni et al., ). Other methods involve learning additional geometric structures, such as cellular sheaves (Bodnar et al., ).

Extending GNNs to directed graphs has been addressed through various methods. GatedGNN (Li et al., ) processed messages from out-neighbors in directed graphs. Some works constructed symmetric matrices from directed adjacency matrices and their transposes to build standard Laplacians (Tong et al., ; Kipf & Welling, ), while others (Ma et al., ; Tong et al., ) developed Laplacians based on random walks and PageRank (Duhan et al., ). MagNet (Zhang et al., ) utilized the magnetic Laplacian to represent directed messages, a technique also applied in adapting transformers to directed graphs (Geisler et al., ). FLODE (Maskey et al., ) employed asymmetrically normalized adjacency matrices within a neural ODE framework. DirGNN (Rossi et al., ) separately processed the messages from in-neighbors and out-neighbors using asymmetrically normalized adjacency matrices, improving node classification on heterophilic graphs. A similar strategy was used in Koke & Cremers (), replacing the adjacency matrices with filters representing Faber polynomials. Recent graph PDE-based models (Eliasof et al., ; Zhao et al., ) introduced an advection term to model directional feature propagation alongside diffusion, assigning edge weights based on computed velocities between nodes, akin to attention coefficients in GAT.

GNNs have been increasingly applied in domains relevant to our work. In single-cell biology, GNNs have been used to predict perturbation responses in gene expression data (Roohani et al., ; Molho et al., ), with datasets compiled in scPerturb (Peidli et al., ). In web traffic analysis, a form of spatiotemporal data on graphs, GNNs often model temporal signals using recurrent neural networks on graphs (Li et al., ; Chen et al., ; Sahili & Awad, ), with datasets and benchmarks provided by PyTorch Geometric Temporal library (Rozemberczki et al., ). In power grids, GNNs have been applied to predict voltage values (Ringsquandl et al., ) and solve optimal power flow problems (Donon et al., ; Böttcher et al., ; Piloto et al., ), with datasets compiled by Lovett et al. ().

## 5 Experiments

### 5.1 Node classification

We first benchmarked our method for node classification task on five graphs with varying properties without simultaneously learning the edge directions. While this is not the main goal of our paper, we include this analysis nonetheless to highlight the advantage of our Laplacian over the magnetic Laplacian as well as the benefit of processing self, in-neighbor aggregated, and out-neighbor aggregated features separately. Importantly, we do not learn edge directions in this case, and hence the phase value is either 0 0 or $\pi/2$ for directed graphs and $\pi/4$ for undirected graphs. For comparison, we include the classical models: GCN (Kipf & Welling, ), SAGE (Hamilton et al., ), GAT (Veličković et al., ); heterophily-specific model, GGCN (Yan et al., ); directionality-aware models: MagNet (Zhang et al., ), FLODE (Maskey et al., ), DirGNN (Rossi et al., ); and a model that learns geometric structure of graph, Sheaf (Bodnar et al., ). We also include MLP to highlight the effect of solely processing the nodes’ self-features without aggregating features across the graph.

As shown in Table 1, CoED demonstrates competitive performance across all five datasets, ranking within the top three in terms of test accuracy. While all models exhibit comparable results on the Cora graph—which is undirected and homophilic—their performances differ significantly on the directed, heterophilic graphs. The classical models developed for undirected graphs particularly struggle on these datasets, with the exception of SAGE, which achieves over 80% accuracy on Texas and Wisconsin. This is because processing only the node’s own features yields good performance, as evidenced by the MLP’s results. In contrast, for the Squirrel and Chameleon datasets, processing directed messages along only one direction is crucial for good performance. Only FLODE, DirGNN, and CoED exhibit strong results on these datasets when configured accordingly. Specifically, for CoED, we introduce the $\alpha$ hyperparameter as in Rossi et al. () to weigh the directional messages post aggregation and make the transformation of self-features optional. The selected hyperparameters modify Equation 4 to $\sigma(\mathbf{m}^{(l)}_{\leftarrow}\mathbf{W}_{\leftarrow}^{(l)}+\mathbf{B}^{
(l)})$ for these datasets.

Our results highlight the advantage of the fuzzy Laplacian over the magnetic Laplacian. The magnetic Laplacian does not process the aggregated messages from out-neighbors and in-neighbors separately. Instead, it combines them into both the real and imaginary components of the aggregated feature vector, thus losing the opportunity to process the two separately. Moreover, during the Laplacian convolution, the directed messages further mix with self-features encoded in the real component. This results in poor performance by MagNet on Squirrel and Chameleon. Sheaf also suffers on these datasets despite expanding the feature dimensions via an object called a stalk, because the sheaf Laplacian is constrained to be symmetric, thereby losing the ability to process directed messages. FLODE and DirGNN, however, perform poorly on Texas and Wisconsin because they focus exclusively on directed messages and do not process the node’s own features separately. Taken together, our benchmarking demonstrates that CoED’s ability to effectively process self-features and separately aggregate in-neighbor and out-neighbor messages using our fuzzy Laplacian enables it to achieve competitive performance across diverse datasets.

**Table 1: Comparison of baseline models and CoED (without edge direction learning) for node classification task across different types of graphs. Top three models are colored by First, Second, Third. The reported numbers are the mean and standard deviation of test accuracies across different splits. The first two rows report the homophily ratios of the graphs and whether they are directed or undirected.**
|  | Texas | Wisconsin | Squirrel | Chameleon | Cora |
| --- | --- | --- | --- | --- | --- |
| Hom. level | 0.11 | 0.21 | 0.22 | 0.23 | 0.81 |
| Undirected | ✗ | ✗ | ✗ | ✗ | ✓ |
| MLP | 80.81$\pm$4.75 | 85.29$\pm$3.31 | 28.77$\pm$1.56 | 46.21$\pm$2.99 | 75.69$\pm$2.00 |
| GCN | 55.14$\pm$5.16 | 51.76$\pm$3.06 | 53.43$\pm$2.01 | 64.82$\pm$2.24 | 86.98$\pm$1.27 |
| SAGE | 82.43$\pm$6.14 | 81.18$\pm$5.56 | 41.61$\pm$0.74 | 58.73$\pm$1.68 | 86.90$\pm$1.04 |
| GAT | 52.16$\pm$6.63 | 49.41$\pm$4.09 | 40.72$\pm$1.55 | 60.26$\pm$2.50 | 86.33$\pm$0.48 |
| GGCN | 84.86$\pm$4.55 | 86.86$\pm$3.29 | 55.17$\pm$1.58 | 71.14$\pm$1.84 | 87.95$\pm$1.05 |
| MagNet | 83.3$\pm$6.1 | 85.7$\pm$3.2 | 39.01$\pm$1.93 | 58.22$\pm$2.87 | 82.63$\pm$1.80 |
| FLODE | 77.57$\pm$5.28 | 80.20$\pm$3.56 | 74.03$\pm$1.58 | 77.98$\pm$1.05 | 86.44$\pm$1.17 |
| Sheaf | 85.95$\pm$5.51 | 89.41$\pm$4.74 | 56.34$\pm$1.32 | 68.68$\pm$1.73 | 87.30$\pm$1.15 |
| DirGNN | 68.38$\pm$ 4.99 | 68.82$\pm$ 5.91 | 75.31$\pm$1.92 | 79.71$\pm$1.26 | 86.27$\pm$ 1.45 |
| CoED | 84.59$\pm$4.53 | 87.84$\pm$3.70 | 75.32$\pm$1.82 | 79.69$\pm$1.35 | 87.02$\pm$1.01 |

### 5.2 Node regression on graph ensemble dataset

Our key contribution is the joint learning of continuous edge directions alongside GNN parameters. This approach is most effective on graph ensemble data, where the graph structure remains fixed but multiple realizations of node features and targets exist. By learning edge directions for all edges without a need to mask parts of the graph, our method optimizes information flow across the entire graph. We empirically demonstrate the benefits of learning continuous edge directions using both synthetic and real-world graph ensemble datasets.

#### 5.2.1 Synthetic datasets

##### Directed flow on triangular lattice.

We begin by applying CoED GNN to a node regression problem constructed on a graph with continuous edge directions, where the target node features are obtained by directionally message-passing the input node features over long distances across the graph. To generate such a graph with continuous edge directions exhibiting long-range order, we created a two-dimensional triangular lattice, assigning each node a position in the 2d plane. We then defined a potential energy function $V$ on this plane, consisting of one peak and one valley (Figure 3(a)). The gradient of $V$ yields a vector field with long-range order, which we used to assign continuous edge directions to the edges of the triangular lattice (Figure 3(b)). Using this graph, we performed the message passing step of Equation 4 iteratively 10 times—using random matrices $\mathbf{W}_{\rightarrow}$, $\mathbf{W}_{\leftarrow}$, and $\mathbf{W}_{\text{self}}$ that were shared across all 10 iterations, starting from the initial node features to obtain the target node values. We repeated this procedure 500 times for different random initial node features and generated an ensemble of input node features and corresponding target node values. During training, we provided all models with the undirected version of the triangular lattice graph (i.e., all $\theta_{ij}=\pi/4$ for CoED). The goal of the learning task is to predict the target node values from the input node features. Additionally, CoED GNN is expected to learn the underlying ground truth continuous edge directions of the graph as part of its training. Further details on data generation are provided in Appendix A.2.1.

##### Gene Regulatory Network (GRN) dynamics.

Gene Regulatory Networks (GRNs) are directed graphs where nodes represent genes and edges represent interactions between pairs of genes (Figure 3(c)). In these networks, when two genes interact, one either activates or suppresses the other. We used Hill functions with randomly chosen parameters to define the dynamics of these gene-gene interactions. We constructed a directed GRN graph with 200 nodes and randomly assigned interactions between them. Starting from random initial expression levels, we solved the system of nonlinear ordinary differential equations representing the GRN dynamics to obtain the steady-state expression levels of all genes. Next, we modeled gene perturbations by setting the expression levels of either one or two genes (the perturbed set of genes) to zero and recomputing the steady-state expression levels for all genes using the same GRN dynamics (Figure 3(d)). We performed this procedure for all single-gene perturbations and a subset of double-gene perturbations, resulting in 1,200 different realizations. Our learning task is to predict the steady-state expression levels of all genes following perturbation (target node values) given the initial steady-state expression levels, with the perturbed genes set to zero (input node features). We provided baseline models with the original graph and CoED with the undirected version of the graph. Further details of the data generation are provided in Appendix A.2.2.

Figure: Figure 3: Synthetic datasets. (a-b) Triangular lattice graph with edge directions derived from the gradient of a 2D potential function $V$ (shown in a), creating long-range flows across the graph. (c-d) Gene regulatory network (GRN) represented as a directed graph where nodes are genes and edges denote interactions. Steady-state gene expression levels are obtained from GRN dynamics, with perturbations simulated by setting the expression levels of specific genes to zero.
Refer to caption: x3.png

##### Results.

As shown in Table 5, CoED achieves the best performance on both synthetic datasets by a considerable margin. For comparison, we include several baseline models: classical models, GCN and GAT; a transformer-based model with positional encoding, GraphGPS (Rampášek et al., ); a directionality-aware model, MagNet; a higher-order model, DRew (Gutteridge et al., ); and a combination of directionality-aware and higher-order model, FLODE. Details of the training setup, hyperparameter search procedure, and selected hyperparameters are provided in Appendix A.2.2.

To identify which aspects of GNNs are particularly effective for learning on graph ensemble datasets, we analyze the baseline models’ results in detail. On the undirected lattice graph, MagNet provides only a slight improvement over GCN, which is expected since MagNet reduces to ChebNet (Defferrard et al., ) on undirected graphs, and GCN is a first-order truncation of ChebNet. However, in the directed GRN experiments, MagNet shows substantial improvement over GCN. We also observe that higher-order GNNs like DRew and FLODE perform competitively on the undirected lattice graph. Notably, FLODE’s instantaneous enhancement of connectivity via fractional powers of the graph Laplacian outperforms DRew’s more gradual incorporation of higher-hop messages. However, both methods struggle on the directed GRN graph. On both synthetic graphs, attention-based models—GAT and GraphGPS—deliver strong performance, coming in just behind CoED. GraphGPS, in particular, seems to benefit from its final global attention step, similar to how FLODE benefits from densifying the graph. We also notice that increasing the dimension of Laplacian positional encoding does not further enhance GraphGPS’s performance. Interestingly, models that learn edge weights via attention mechanisms outperform MagNet on the directed GRN graph. This is likely because MagNet’s unitary evolution of complex-valued features does not resemble the actual feature propagation (i.e., the GRN dynamics), in addition to the shortcomings highlighted in the previous section.

We then investigated whether CoED can recover the true continuous edge directions of the triangular lattice graph, given that the feature propagation steps during data generation closely match the message-passing operation of CoED. As shown in Figure 5, CoED correctly learns the true directions. Lastly, since both synthetic datasets are generated by propagating input features over multiple hops, we investigated how performance scales with model depth by training CoED and the second-best model with up to 10 layers. Figure 6 demonstrates that CoED continues to improve as depth increases, while the performance of the other models plateau at a shallower depth.

Figure: Figure 4: Comparison of different models on the synthetic datasets. Values are test losses reported with a common factor of $10^{-3}$ in both columns.
Refer to caption: extracted/5934640/figs/fig4.png

Figure: Figure 6: Model performance as a function of depth.
Refer to caption: x4.png

#### 5.2.2 Real Datasets

##### Single-cell Perturb-seq.

Perturb-seq (Dixit et al., ) is a well-established experimental technique in single-cell biology that inspired the synthetic GRN experiment described earlier. In Perturb-seq experiments, one or more genes in a cell are knocked out resulting in zero expresion—as in our synthetic GRN dataset. The resulting changes in the expression levels of all other genes are then measured to elucidate gene-gene interactions. For our study, we used the Replogle-gwps dataset (Replogle et al., ; Peidli et al., ), which includes 9,867 distinct single-gene perturbations, along with control measurements from cells without any perturbation to establish baseline gene expression levels. The learning task is again predicting the expression levels of all genes following perturbation given the initial steady-state with the expression levels of the perturbed genes set to zero. Since there is no ground truth gene regulatory network (GRN) available for this dataset, we constructed an undirected $k$-nearest neighbors graph to connect genes with highly correlated expression patterns. All models are trained using this heuristic graph. Details of the data processing procedure are provided in Appendix A.3.1.

##### Wikipedia web traffic.

We also modeled the traffic flow between Wikipedia pages using the WikiMath dataset, which is classified as a “static graph with temporal signals” in the PyTorch Geometric Temporal library (Rozemberczki et al., ). In this dataset, each node corresponds to an article page on a popular mathematics topic, and each directed edge represents a link from one page to another. The node features are the daily visit counts of all pages over a period of 731 consecutive days. The learning task is node regression: predict the next day’s visit counts across all pages given today’s visit counts. We trained the baseline models using the true, directed graph, while CoED was trained starting from the undirected version of the graph. Additional details are provided in Appendix A.3.2.

##### Power grid.

We applied CoED to the optimal power flow (OPF) problem using the OPFData (Lovett et al., ) from the PyTorch Geometric library. In this dataset, a power grid is represented as a directed graph with nodes corresponding to buses (connection points for generators and loads) and edges representing transformers and AC lines. Input features are the operating values of all components under specific load conditions, and the targets are the corresponding AC-OPF solution values at the generator nodes. To compare different models, we used a consistent architecture across components but substituted different model layers for message passing. For CoED, they were again converted to undirected edges. Additional details are provided in Appendix A.3.3.

##### Results.

**Table 2: Comparison of different methods on real graph ensemble datasets. Values are test losses reported with common factors of $10^{1},10^{-1},10^{-3}$ for Perturb-seq, web traffic, and power grid columns, respectively.**
|  | Perturb-seq | Web traffic | Power grid |
| --- | --- | --- | --- |
| GCN | 4.13$\pm$0.08 | 7.07$\pm$0.03 | 28.56$\pm$6.08 |
| MagNet | 4.11$\pm$0.01 | 6.94$\pm$0.02 | 18.05$\pm$2.77 |
| GAT | 3.85$\pm$0.03 | 6.00$\pm$0.03 | 13.57$\pm$1.73 |
| DirGCN | 5.46$\pm$0.26 | 6.72$\pm$0.04 | 6.15$\pm$0.84 |
| DirGAT | 3.98$\pm$0.07 | 6.55$\pm$0.04 | 3.28$\pm$0.17 |
| CoED | 3.56$\pm$0.03 | 5.76$\pm$0.05 | 2.91$\pm$0.11 |

Table 2 reports the test performances of all baseline models and CoED on the three datasets. The baselines include GAT, MagNet, DirGCN, and DirGAT. We focus on these models because attention-based approaches showed competitive performance on the synthetic datasets, and MagNet, which accounts for edge directions, performed well on the directed GRN dataset. Details of the training setup, hyperparameter search procedure, and selected hyperparameters are provided in the Appendix.

We observe that CoED achieves the best performance across all three datasets. On the Perturb-seq dataset with an undirected graph, MagNet performs similarly to GCN, while DirGCN struggles. We attribute DirGCN’s poor performance to clashing learnable parameters: it uses two distinct weight matrices, $\mathbf{W}_{\leftarrow}$ and $\mathbf{W}_{\rightarrow}$, applied to identical in- and out-neighbor aggregated messages in the case of undirected graphs. For a propagation path of $L$ hops, this results in $2^{L}$ feature transformations, comprised of different combinations of the two weight matrices, which together reduce the model’s ability to efficiently learn the optimal weight matrices. In contrast, DirGAT’s attention mechanisms break the symmetry of the undirected edges, leading to improved performance. CoED naturally addresses this issue in undirected graphs by learning the edge directions. On the web traffic and power grid datasets, which have directed graphs, we observe a similar trend. MagNet outperforms GCN due to its ability to process directed messages. However, GAT delivers better performance than MagNet, likely because its attention mechanism effectively captures important features. Since directed graphs create distinct feature propagation paths, DirGCN achieves substantial performance gains. DirGAT further improves upon DirGCN by leveraging additional edge weight learning through an attention mechanism. CoED surpasses all these models, demonstrating the effectiveness of learning continuous edge directions.

## 6 Conclusion

We have introduced the Continuous Edge Direction (CoED) GNN, which assigns fuzzy, continuous directions to the edges of a graph and employs a novel complex-valued Laplacian to transform information propagation on graphs from diffusion to directional flow. Our theoretical analysis shows that CoED GNN is more expressive than existing Laplacian-based methods and matches the expressiveness of an extended Weisfeiler-Leman test for directed graphs with fuzzy edges. Through extensive experiments on both synthetic and real-world graph ensemble datasets—including gene regulatory networks, web traffic, and power grids—we demonstrated that learning continuous edge directions significantly improves performance over existing GNN models. These results highlight CoED GNN’s effectiveness in enabling long-range information flow, offering a powerful framework for learning on graphs.

## Appendix A Experimental details

### A.1 Node classification

All the results reported in Table 1 are taken from (Bodnar et al., ) except for MagNet, FLODE, and DirGNN. Magnet’s results on Texas and Wisconsin are taken from the original paper, while those on Squirrel and Chameleon are taken from (Rossi et al., ). FLODE’s results on Squirrel, Chameleon, and Cora, and DirGNN’s results on Squirrel and Chameleon are taken from their respective papers. We trained these models and CoED for the remaining results and describe the training procedures below. All five datasets were downloaded using PyTorch Geometric library (Fey & Lenssen, ) with split=‘geom-gcn’ argument to use the 10 fixed 48%/32%/20% training/validation/test splits provided by (Pei et al., ), which are the splits used in all the referenced results.

##### Training.

We evaluated the validation accuracy at each epoch, incrementing a counter if the value did not improve and resetting it to 0 when a new best validation accuracy was achieved. Training was early-stopped when the counter reached a patience of 200. Unless otherwise mentioned, we used the default hyperparameter settings of the respective models. For DirGNN, we chose GCN as a convolution layer since it showed competitive performance across various datasets considered in the original paper. We used the ReLU activation function and the ADAM optimizer in all experiments. Across all models, we searched over the following hyperparmeters: hidden dimension $\in\{32,64,128,256\}$, learning rate $\in\{\text{5e-4},\text{1e-3},\text{2e-3},\text{5e-3},\text{1e-2},\text{2e-2}\}$, and weight decay $\in\{0,\text{1e-4},\text{5e-4},\text{1e-3}\}$. With the exception of FLODE, we set the dropout rate to 0.5 and tuned the number of layers $\in\{2,3,4,5\}$. We additionally searched over model-specific hyperparameters: the weight between in-/out-neighbor aggregated messages $\alpha\in\{0,0.5,1\}$, jumping knowledge (jk) $\in\{\texttt{None},\text{`cat'},\text{`max'}\}$, and layer-wise feature normalization (norm) $\in\{\texttt{True},\texttt{False}\}$ for DirGNN and CoED; additionally, self-feature transform $\in\{\texttt{True},\texttt{False}\}$, self-loop value $\in\{0,1\}$ for CoED; the order of Chebyshev polynomial $K\in\{1,2\}$, the global directionality $q\in\{0,0.05,0.1,0.15,0.2,0.25\}$, and self-loop value $\in\{0,1\}$ for MagNet; and the number of layers $\in\{1,2,3\}$, self-loop value $\in\{0,1\}$, and dropout rate for encoder and decoder $\in\{0,0.5\}$ for FLODE. Note that the number of layers refers to the number of forward Euler steps in FLODE. We set the number of layers in both encoder and decoder MLPs to 1 to focus on the effect of the fractional Laplacian, and solved the heat equation with a minus sign, which was the setup employed for Squirrel and Chameleon experiments in the original paper. If the self-loop value is 1, the self-feature is combined with neighbors’ features in the AGGREGATE function. We report in Table 1 the mean accuracy and standard deviation over the 10 test splits using the best hyperparameters presented below. All experiments were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB of memory and it took roughly 5 days of training to produce the results.

**Table A.1: Hyperparameters selected for node classification experiments**
|  | # layers | # hidden | lr | wd | dropout | self-loop | K | q | $\alpha$ | norm | jk | self-feature |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| MagNet (Cora) | 1 | 32 | 1e-2 | 0 | 0.5 | 1 | 2 | 0 | - | - | - | - |
| FLODE (Texas) | 1 | 128 | 1e-2 | 1e-3 | 0 | 0 | - | - | - | - | - | - |
| FLODE (Wisconsin) | 1 | 128 | 5e-3 | 0 | 0.5 | 0 | - | - | - | - | - | - |
| DirGNN (Cora) | 2 | 64 | 5e-3 | 5e-4 | 0.5 | - | - | - | 0 | False | None | - |
| DirGNN (Texas) | 4 | 64 | 5e-3 | 1e-3 | 0.5 | - | - | - | 0.5 | False | None | - |
| DirGNN (Wisconsin) | 2 | 128 | 5e-3 | 1e-4 | 0.5 | - | - | - | 0.5 | False | None | - |
| CoED (Texas) | 2 | 64 | 2e-2 | 5e-4 | 0.5 | 0 | - | - | 0.5 | False | None | True |
| CoED (Wisconsin) | 2 | 128 | 2e-2 | 1e-3 | 0.5 | 0 | - | - | 0.5 | False | None | True |
| CoED (Squirrel) | 3 | 128 | 1e-2 | 0 | 0 | 0 | - | - | 0 | True | max | False |
| CoED (Chameleon) | 4 | 128 | 5e-3 | 0 | 0 | 0 | - | - | 0 | True | cat | False |
| CoED (Cora) | 2 | 128 | 5e-4 | 1e-4 | 0.5 | 1 | - | - | 0 | False | None | False |

### A.2 Graph ensemble experiment with synthetic datasets

#### A.2.1 Data generation and training setup for the directed flow triangular lattice graph

##### Data generation.

To obtain the triangular lattice graph described in the main text, we first designed a potential function $V$ on $[-2,2]^{2}$ plane with a peak (source) and a valley (sink). We used quadratic potentials, located at $\mu_{1}=(-1,1)$ and $\mu_{2}=(1,-1)$ with stiffness matrices,

$$ $\mathbf{K}_{1}=\mathbf{K}_{2}=\begin{pmatrix}1&0\\ 0&1\end{pmatrix}$ $$

With magnitudes $a_{1}=1$ and $a_{2}=-1$, the potential function $V(\mathbf{x})$ is parameterized as,

$$ $V(\mathbf{x})=a_{1}(\mathbf{x}-\mu_{1})^{\top}\mathbf{K}_{1}(\mathbf{x}-\mu_{1 })+a_{2}(\mathbf{x}-\mu_{2})^{\top}\mathbf{K}_{2}(\mathbf{x}-\mu_{2})$ $$

We then generated a triangular lattice on the 2d plane and considered this lattice as a graph $\mathcal{G}_{\text{lattice}}$, where the vertices of the lattice serve as the nodes of the graph, and the edges of the lattice form the edges of the graph. In this way, each node $v_{i}\in\mathcal{V}_{\text{lattice}}$ has an associated spatial position $\mathbf{x}_{i}$ on the 2d plane. We computed $\Delta V_{ij}=V(\mathbf{x}_{j})-V(\mathbf{x}_{i})$ for all $(v_{i},v_{j})\in\mathcal{E}_{\text{lattice}}$. All $\Delta V_{ij}$ values were shifted and scaled to the range $[0,\pi/2]$ to obtain the $\theta_{ij}$, which is an approximate version of the gradient direction of the potential. In the resulting lattice graph, an edge points towards the node with the lower potential energy.

We then assigned to each node a 10-dimensional random feature vector sampled independently from the standard multivariate normal distribution, and normalized them to have a unit-norm, and repeated the process 500 times to generate an ensemble of node features. To generate corresponding target values, we propagated features using Equation 4 with message-passing matrices $\mathbf{P}_{\rightarrow}$ and $\mathbf{P}_{\leftarrow}$ computed from $\Theta$ of the lattice graph as described in Equation 2 and the entries of the $10\times 10$ weight matrices $\mathbf{W}_{\rightarrow}$, $\mathbf{W}_{\leftarrow}$, and $\mathbf{W}_{\text{self}}$ sampled independently from the standard normal distribution and shared across all 10 iterations. Instead of applying an activation function, we normalized the features $\mathbf{m}_{\text{self}}+\mathbf{m}_{\rightarrow}+\mathbf{m}_{\leftarrow}$ to have unit norm. We used the features after 10 iterations of message passing as the target values. We divided these 500 instances of feature-target pairs using 60%/20%/20% random training/validation/test split.

##### Training.

We used a batch size of 16 for training with random shuffling at each epoch and a full batch for both validation and testing. We evaluated the validation MSE at each epoch, incrementing a counter if the value did not improve and resetting it to 0 when a new best validation MSE was achieved. Training was early-stopped when the counter reached a patience of 20. We used neither dropout nor weight decay, as we aim to learn an exact mapping from node features to target values for regression, as opposed to a noise-robust node embedding for classification. We used the ReLU activation function in all models, except for ELU in GAT, and used ADAM optimizer for all experiments. Across all models, we searched over the following hyperparmeters: the number of layers $\in\{2,3,4\}$, hidden dimension $\in\{16,32,64\}$, learning rate $\in\{\text{1e-3},\text{2e-3},\text{5e-3},\text{1e-2}\}$. We additionally grid-searched over model-specific hyperparameters: the number of attention heads $\in\{1,4,8\}$ and skip connection (sc) $\in\{\texttt{True},\texttt{False}\}$ for GAT; attention type $\in\{\text{`multihead'},\text{`performer'}\}$, attention heads $\in\{1,4,8\}$, encoding type $\in\{\text{`eigenvector'},\text{`electrostatic'}\}$, the dimension of eigenvector encoding $\in\{2,5,10,20\}$ and self-loop value $\in\{0,1\}$ for GraphGPS; the order of Chebyshev polynomial $K\in\{1,2\}$ and self-loop value $\in\{0,1\}$ for MagNet; multi-hop aggregation mechanism $\in\{\text{`sum'},\text{`weight'}\}$ for DRew; the number of encoder MLP layers $\in\{1,2,3\}$, the number of decoder MLP layers $\in\{1,2,3\}$ and self-loop value $\in\{0,1\}$ for FLODE; self-feature transform $\in\{\texttt{True},\texttt{False}\}$, learning rate for $\Theta\in\{\text{1e-3},\text{5e-3},\text{1e-2}\}$, and layer-wise (lw) $\Theta$ learning $\in\{\texttt{True},\texttt{False}\}$ for CoED. For GraphGPS, we used GINE as a convolution layer and provided 1 as an edge attribute. Computing structural encoding via random walk resulted in an all-zero vector since the degree of node is 3 across all nodes except at the boundary in our lattice graph. We thus opted to use the electrostatic function encoding provided in the original paper as an alternative to structural encoding. For MagNet, we optimized $q$ along with the model parameters during training. We supplied the undirected version of the lattice graph by setting all $\theta_{ij}$ values to $\pi/4$ and used self-feature transform with self-loop value set to 0. We did not use layer-wise feature normalization. Table 5 reports the mean accuracy and standard deviation on the test data from the top 5 out of 7 training runs with different initializations using the best hyperparameters shown below. All experiments were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB of memory and it took about 3 days of training time to generate the results.

**Table A.2: Hyperparameters selected for node regression on the synthetic lattice graph.**
| Model | # layers | # hidden | lr | self-loop | sc | # attn. heads | attn. type | enc. type | K | aggr. | # enc. layers | # dec. layer | self-feature | lr $\Theta$ | lw $\Theta$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| GCN | 2 | 16 | 1e-3 | - | - | - | - | - | - | - | - | - | - | - | - |
| GAT | 4 | 32 | 1e-3 | - | True | 4 | - | - | - | - | - | - | - | - | - |
| GraphGPS | 4 | 64 | 1e-3 | 0 | - | 4 | multihead | eigenvector (dim=5) | - | - | - | - | - | - | - |
| MagNet | 4 | 64 | 1e-3 | 1 | - | - | - | - | 2 | - | - | - | - | - | - |
| DRew | 3 | 64 | 5e-3 | - | - | - | - | - | - | weight | - | - | - | - | - |
| FLODE | 4 | 64 | 1e-2 | 0 | - | - | - | - | - | - | 1 | 3 | - | - | - |
| CoED | 4 | 64 | 1e-3 | - | - | - | - | - | - | - | - | - | True | 1e-3 | False |

#### A.2.2 Data generation and Training setup for the GRN dynamics experiment

##### Data generation.

We prepared a directed adjacency matrix $\mathbf{A}$ for a graph with $200$ nodes by sampling each entry of the matrix independently from a Bernoulli distribution with success probability 0.03. We interpreted an edge $A_{ij}$ as indicating that the gene represented by the node $v_{j}$ is regulating the gene represented by the node $v_{i}$. We then randomly chose half of the edges as activating edges and the other half as suppressing edges. The scalar feature value of each node is the expression level of a gene, measured as concentration $c_{i}$. In order to simulate the gene regulatory network dynamics where genes are either up- or down- regulating connected genes, we sampled the magnitudes of activation $\gamma^{\text{act}}_{ij}$ and suppression $\gamma^{\text{sup}}_{ij}$ from a uniform distribution with the support $[0.5,1.5]$. We additionally sampled half-saturation constants $K_{ij}$, which control how quickly $c_{i}$ changes in response to $c_{j}$, from a uniform distribution with support $[0.25,0.75]$. Lastly, we sampled the initial concentrations independently from a uniform distribution with support $[0.1,10]$ for each gene, and ran the GRN dynamics as described by,

$$ $\frac{dc_{i}}{dt}=\sum_{j\in\mathcal{N}(i)}\big{(}\gamma^{\text{act}}_{ij}F^{ \text{act}}(c_{j},K_{ij})+\gamma^{\text{sup}}_{ij}F^{\text{sup}}(c_{j},K_{ij}) \big{)}-c_{i}$ (5) $$

for 250 time steps with $dt=0.05$ to reach a steady state for $\mathbf{c}$. The summation in the above equation is over all genes that either activate or repress gene $i$. $F^{\text{act}}(c_{j},K_{ij})=c_{j}^{2}/(K_{ij}^{2}+c_{j}^{2})$ and $F^{\text{sup}}(c_{j},K_{ij})=K_{ij}^{2}/(K_{ij}^{2}+c_{j}^{2})$, which are the Hill functions defining up- and down- regulations of gene $i$ by gene $j$, respectively. From this steady state, we mimicked the gene knockout experiment in biology by setting the concentration values of a chosen set of genes to zero and running the GRN dynamics for additional 100 time steps, by which time genes reached new steady-state values. We performed a single-gene knockout for all 200 genes and a double-gene knockout for 1000 randomly selected pairs of genes. We defined node features as the original steady state with the values of knockout genes set to zero and the corresponding target values as the new steady-state values reached from this state. This procedure generates an ensemble of 1200 feature-target pairs for each node for the synthetic GRN graph. We used all 200 single gene knockout results as training data, and randomly selected 200 and 800 double gene knockout results for validation and testing, respectively.

##### Training.

For each knockout result, we used all nodes for regression except those corresponding to the knocked out genes. We used a batch size of 8 for training with random shuffling at each epoch, and a full batch for both validation and testing. We evaluated the validation loss at every epoch, and implemented the same counting scheme as in the directed flow experiment to early-stop the training with a patience of 50. We searched over the number of layers $\in\{2,3,4,5\}$, hidden dimension $\in\{16,32\}$, and learning rate $\in\{\text{5e-4},\text{1e-3},\text{2e-3},\text{5e-3}\}$, and otherwise conducted the same hyperparameter search as described in the lattice experiment, using the same training setup. Table 5 reports the mean accuracy and standard deviation on the test data from the top 5 out of 7 training runs with different initializations using the best hyperparameters shown below. All experiments were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB of memory and it took about 3 days of training time to generate the results.

**Table A.3: Hyperparameters selected for node regression on the synthetic GRN graph.**
| Model | # layers | # hidden | lr | self-loop | sc | # attn. heads | attn. type | enc. type | K | aggr. | # enc. layers | # dec. layer | self-feature | lr $\Theta$ | lw $\Theta$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| GCN | 3 | 32 | 5e-4 | - | - | - | - | - | - | - | - | - | - | - | - |
| GAT | 5 | 32 | 2e-3 | - | True | 8 | - | - | - | - | - | - | - | - | - |
| GraphGPS | 5 | 32 | 5e-3 | 0 | - | 4 | multihead | eigenvector (dim=10) | - | - | - | - | - | - | - |
| MagNet | 5 | 32 | 5e-3 | 1 | - | - | - | - | 2 | - | - | - | - | - | - |
| DRew | 3 | 32 | 5e-4 | - | - | - | - | - | - | weight | - | - | - | - | - |
| FLODE | 5 | 32 | 5e-3 | 1 | - | - | - | - | - | - | 1 | 1 | - | - | - |
| CoED | 5 | 32 | 1e-3 | - | - | - | - | - | - | - | - | - | True | 1e-2 | True |

### A.3 Graph Ensemble experiment with real datasets

#### A.3.1 Preprocessing, Data generation, and training setup for single-cell perturbation experiments

##### Preprocessing.

We downloaded Replogle-gwps dataset from scPerturb database (Peidli et al., ) and followed the standard single-cell preprocessing routine using Scanpy software (Wolf et al., ), selecting for the top 2000 most variable genes. This involved running the following four functions: filter_cells function with min_counts=20000 argument, normalize_per_cell function, filter_genes function with min_cells=50 argument, highly_variable_genes function with n_top_genes=2000, flavor=‘seurat_v3’, and layer=‘counts’ arguments. Afterwards, we discarded genes that were not part of the top 2000 most variable genes. These 2000 genes define nodes. We did not log transform the expression values and used the normalized expression values obtained from these preprocessing steps for all downstream tasks. Out of 9867 genes that were perturbed in the original dataset, 958 of them were among the top 2000 most variable genes. Note that perturbed genes should have near zero expression values since the type of perturbation in the original experiment was gene knockout via a technique called CRISPRi. Therefore, we disregarded perturbed genes if the perturbations did not result in more than 50% of cells with zero expression values for each respective gene. This preprocessing step identifies 824 effective gene perturbations. Note that there are multiple measurements per perturbation.

##### Data generation for node regression.

In perturbation experiments, the expression values are measured ‘post-perturbation’ (equivalent to $c_{i}+\Delta c_{i}$ in Figure 3(d) of the main text). Consequently, we do not have access to their ‘pre-interaction’ expression levels (equivalent to $c_{i}$ with a perturbed gene’s value set to 0). To pair each post-perturbation expression values to its putative pre-interaction state, we randomly sampled a control measurement and set the expression value of the perturbed gene to 0. These pairs of pre-interaction and post-perturbation expressions define the ensemble of features and targets. Since the number of measurements vary per perturbation, we standardized the dataset diversity by downsampling the feature target pairs to 2 per perturbation. We split each dataset based on perturbations, grouping all cells subject to the same perturbation together. We performed 60/10/30 training/validation/test splits and computed $k$-nearest neighbors gene-gene graph with $k=3$ using the training split. We checked that this graph roughly corresponds to creating an edge between genes whose Pearson correlation coefficient is higher than 0.5. We also confirmed that this graph represents a single connected component.

##### Training.

As in the GRN example, we used all nodes (i.e., genes) for regression except for those corresponding to a knocked-out gene. We used a batch size of 16 for training with random shuffling at each epoch, and a full batch for both validation and testing. We evaluated the validation loss at every epoch, and implemented the same counting scheme as in the synthetic dataset experiments to early-stop the training with a patience of 30. Since the models’ performances generally improved as their depths and hidden dimension increased, we used $4$ layers and a hidden dimension of 32 across all models (to streamline the comparison). We searched over learning rate $\in\{\text{1e-3},\text{5e-3},\text{1e-2}\}$ for all models, the order of Chebyshev polynomial $K\in\{1,2\}$ for MagNet; the number of attention heads $\in\{1,4,8\}$ for both GAT and DirGAT; addtionally, skip connection (sc) $\in\{\texttt{True},\texttt{False}\}$ for GAT; self-feature transform $\in\{\texttt{True},\texttt{False}\}$, learning rate for $\Theta$ $\in\{\text{5e-4},\text{1e-3},\text{5e-3},\text{1e-2}\}$ and layer-wise $\Theta$ learning $\in\{\texttt{True},\texttt{False}\}$ for CoED. We learned $q$ in MagNet. Table 2 reports the mean accuracy and standard deviation on the test data from the top 5 out of 7 training runs with different initializations using the best hyperparameters shown below. All experiments were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB memory and approximately one day of training time was spent for generating the results.

**Table A.4: Hyperparameters selected for node regression on the Perturb-seq data.**
| Model | lr | K | sc | # attn. heads | self-feature | lr $\Theta$ | lw $\Theta$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| GCN | 5e-3 | - |  | - | - | - | - |
| MagNet | 1e-3 | 2 |  | - | - | - | - |
| GAT | 1e-3 | - | True | 4 | - | - | - |
| DirGCN | 5e-3 | - |  | - | - | - | - |
| DirGAT | 1e-3 | - |  | 1 | - | - | - |
| CoED | 1e-3 | - |  | - | False | 5e-4 | False |

#### A.3.2 Training setup for web traffic experiments

We downloaded WikiMath dataset from PyTorch Geometric Temporal library (Rozemberczki et al., ) with the default time lag value of 8. We followed the same temporal split from the paper where 90% of the snapshots were used for training and 10% of the forecasting horizons were used for testing.

##### Training.

The input features and target values are shaped as $N\times 8$ and $N\times 1$, respectively. Since we consider a pair of the two consecutive snapshots as input and the corresponding target, we disregarded the first 7 snapshots (i.e., feature dimensions) and used only the node values (i.e., visit counts of Wikipedia articles) of the last snapshot for prediction for testing. For training, we utilized all 8 snapshots, generating 8 predictions. MSE loss for the predictions from the first 7 snapshots were each evaluated with the next 7 snapshots as target values. This procedure is similar to the ‘incremental training mode’ used to train time-series based models (Rozemberczki et al., ; Eliasof et al., ; Guan et al., ). We used a full batch for both training and testing. We evaluated test loss at each step, and implemented the same counting scheme used above to early-stop training with a patience value of 50. We used 2 layers with a hidden dimension of 16 across all models. We searched over learning rate $\in\{\text{1e-3},\text{5e-3},\text{1e-2},\text{2e-2}\}$ for all models and also over the same model-specific hyperparameters discussed in the Perturb-seq experiments except the number of attention heads, which was limited to 2 due to memory constraints. We learned $q$ for MagNet. Table 2 reports the mean accuracy and standard deviation on the test data from the top 5 out of 7 training runs with different initializations using the best hyperparameters shown below. All experiments were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB memory and approximately 16 hours of training time was spent for generating the results.

**Table A.5: Hyperparameters selected for node regression on the WikiMath data.**
| Model | lr | K | sc | # attn. heads | self-feature | lr $\Theta$ | lw $\Theta$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| GCN | 1e-2 | - |  | - | - | - | - |
| MagNet | 5e-3 | 1 |  | - | - | - | - |
| GAT | 1e-2 | - | True | 2 | - | - | - |
| DirGCN | 2e-2 | - |  | - | - | - | - |
| DirGAT | 5e-2 | - |  | 1 | - | - | - |
| CoED | 5e-3 | - |  | - | True | 1e-2 | False |

#### A.3.3 Training setup for power grid experiments

We downloaded a power grid graph with 2000 nodes from PyTorch Geometric library, compiled via OPFData (Lovett et al., ). We selected the ‘fulltop’ topology option to obtain the information of the entire graph, opposed to ‘N-1’ perturbation option which masks parts of the graphs. We randomly sampled 300 graphs and split them into 200 training set, 50 validation set, and 50 test set. Power grids are heterogeneous graphs. Refer to Figure 1 of Piloto et al. () for detailed overview of the different components. Importantly, generators, loads, and shunts (‘subnodes’) are all connected to buses (‘nodes’) via edges, and buses are connected to one another via two types of edges, transformers and AC lines. Thus, the load profile across the graph informs the generator subnodes via the bus-to-bus edges. We describe the architectural design choices that we made to accomplish this task while facilitating effective model comparison.

##### Model.

The primary goal is to ensure information flow into those buses connected to generators. We thus substituted different model layers to process messages over the two types of bus-to-bus edges while keeping the rest of the architecture unchanged. The processing steps are:

- 1.
Transform node/subnode features and edge features using an MLP.
- 2.
For each bus node, integrate the features of its subnodes into its own features via GraphConv (i.e., $\mathbf{W}_{1}\mathbf{x}_{i}+\mathbf{W}_{2}\sum_{j\in\mathcal{N}(i)}e_{j,i}
\cdot\mathbf{x}_{j}$
).
- 3.
Incorporate edge features into node features using GINE.
- 4.
Iterate message-passing among bus nodes
- 5.
Decode the features of bus nodes into generator operating point values.

We varied the message-passing mechanism in step 4 by applying the different Aggregate and Update functions of each of the models that we analyzed.

##### Training.

Following the example training routine outlined in Lovett et al. (), we trained models to predict generator active and reactive power outputs and evaluated MSE loss against those values in the AC-OPF solutions. We did not incorporate AC-OPF constraints, as the focus of the experiments was to compare the message-passing capabilities of CoED with other models. We refer to Böttcher et al. (); Piloto et al. () for predicting AC-OPF solutions that satisfy constraints. We used a batch size of 16 during training with random shuffling applied at each epoch. We evaluated the validation loss at every epoch and early-stopped the training with the same counting scheme with a patience of 50. We iterated the step 4 above 3 times (i.e., 3 layers) and used 32 hidden dimension. We searched over learning rate $\in\{\text{5e-4},\text{1e-3},\text{2e-3},\text{5e-3}\}$, and otherwise the same set of hyperparameters considered in the Perturb-seq experiments. We learned $q$ for MagNet. Table 2 reports the mean accuracy and standard deviation on the test data from the top 5 out of 7 training runs with different initializations using the best hyperparameters shown below. All training were performed on two NVIDIA RTX 6000 Ada Generation GPUs with 48GB memory and approximately one day of training time was spent for generating the results.

**Table A.6: Hyperparameters selected for node regression on the AC-OPF data.**
| Model | lr | K | sc | # attn. heads | self-feature | lr $\Theta$ | lw $\Theta$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| GCN | 5e-3 | - |  | - | - | - | - |
| MagNet | 1e-3 | 2 |  | - | - | - | - |
| GAT | 2e-3 | - | False | 4 | - | - | - |
| DirGCN | 1e-3 | - |  | - | - | - | - |
| DirGAT | 1e-3 | - |  | 8 | - | - | - |
| CoED | 5e-4 | - |  | - | False | 5e-3 | True |

## Appendix B Time and Space Complexity

Let $V$ and $E$ denote the numbers of nodes and edges, and let $H$ denote hidden dimension, which we assume stays constant for two consecutive layers. At minimum, the message-passing mechanism outlined in 2 results in the time complexity that scales as $\mathcal{O}(EH+VH^{2})$, where the first term corresponds to aggregating features with dimension $H$ and the second term stems from the matrix-matrix multiplication of node features and a weight matrix. The space complexity is $\mathcal{O}(E+H^{2})$ due to the edges and the weight matrix.

With the sparse form of the phase matrix $\Theta$, CoED incurs an additional $\mathcal{O}(E)$ term both in the time complexity from computing fuzzy propagation matrices $\mathbf{P}_{\leftarrow/\rightarrow}$ and in the space complexity due to storing $\Theta$. Unless layer-wise $\Theta$ learning is employed, this computation happens once and thus only adds minimal overhead. The layer-wise $\Theta$ learning adds an $\mathcal{O}(EL)$ term to the total time and space complexities over $L$ layers. For comparison, however, we note that computing $S$-head attentions incurs an $\mathcal{O}(EH)$ term (or even $\mathcal{O}(EHS)$ term if feature dimension is not divided by $S$) in the time complexity and an $\mathcal{O}(ES)$ term in the space complexity per layer.

## Appendix C Positional encoding using the fuzzy Laplacian

Graph Laplacians can be used to assign a positional encoding to each node of a graph based on the connectivity patterns of the nodes of the graph. Using the fuzzy Laplacian, we can extend positional encoding to include variations in directions of the edges surrounding a node in addition to the connectivity pattern of the graph. To demonstrate the utility of the eigenvectors of the fuzzy Laplacian for positional encoding, we visualize the eigenvectors of the Laplacian computed from the triangular lattice graph (which has a trivial connectivity pattern as shown above by the random walk structural encoding being trivially zero) supplemented with two different sets of edge directions. In the first case, we obtained edge directions from the gradient of source-sink potential function described in A.2.1. These edge directions are visualized in Figure 3(b) of the main text. In the second case, we obtained the edge directions for the same triangular lattice graph, from the following solenoidal vector field,

$$ $F(\mathbf{x})=\big{(}\sin(\pi x)\cos(\pi y),-\cos(\pi x)\sin(\pi y)\big{)}$ $$

which does not have sources or sinks but instead features a cyclic flow. To assign edge directions using the solenoidal vector field, we computed $\theta_{ij}$ as the angle between the unit vector pointing from $v_{i}$ to $v_{j}$ and the vector evaluated at the midpoint of an edge, i.e., $F\big{(}(\mathbf{x}_{2}-\mathbf{x}_{1})/2\big{)}$. All $\theta_{ij}$ were scaled to range from 0 0 to $\pi/2$. As shown in Figure C.1(a-e), the real part of the eigenvector distinguishes the peak and valley from the region in between, whereas the phase of the eigenvector distinguish the peak from the valley. Similarly in Figure C.1(f-j), the real part of the eigenvector highlights the regions adjacent to the cyclic flows. The magnitude of the eigenvector picks out the centers of the four solenoids in the middle. Taken together, the eigenvectors of the fuzzy Laplacian contain positional information based on the directions of the edges as demonstrated here for the same graph but with different edge directions.

Figure: Figure C.1: Visualization of the eigenvector, corresponding to the eigenvalue with largest magnitude, of the fuzzy Laplacian of the triangular lattice whose edge directions are taken from the source-sink potential function described in the main text (top row) and the solenoidal vector field described in this section (bottom row). The original vector fields are shown in the left-most figures. The real and imaginary components, as well as the magnitude and phase, of the eigenvector encode positional information at each node about the direction of the edges surrounding that node.
Refer to caption: x5.png

## Appendix D Mathematical properties of the fuzzy Laplacian

We propose a new graph Laplacian matrix for directed graphs which generalizes to the case of directed graphs with fuzzy edges, where an edge connecting two nodes $A$ and $B$ can take on any intermediate value between the two extremes of pointing from $A$ to $B$ and pointing from $B$ to $A$. We will show that our Laplacian exhibits two useful properties: 1. The eigenvectors of our Laplacian matrix are orthogonal and therefore can be used
as the positional encodings of the nodes of the graph. 2. For any node of the graph, our Laplacian aggregates information from neighbors that send information to the node separately from the neighbors that receive information from the node and is therefore as expressive as a weak form of the Weisfeiler-Leman (WL) graph isomorphism test for directed graphs with fuzzy edges that we define below.

### D.1 Directed graphs with fuzzy edges

We define a directed graph with fuzzy edges as follows. Our definition builds on the standard definition of a graph.

A graph $\mathcal{G}$ is an ordered pair $\mathcal{G}:=(\mathcal{V},\mathcal{E})$ comprising a set $\mathcal{V}$ of vertices or nodes together with a set $\mathcal{E}$ of edges. Each edge is a 2-element subset of $\mathcal{V}$.

- •
$\mathcal{V}$: A finite, non-empty set of vertices
$\mathcal{V}=\{v_{1},v_{2},\ldots,v_{n}\}$
- •
$\mathcal{E}$: A set of edges, each linking two vertices in $\mathcal{V}$
$\mathcal{E}=\{(v_{i},v_{j})\mid v_{i},v_{j}\in\mathcal{V}\}$

To incorporate fuzzy directions to the edges, we define a new attribute for each edge.

- •
$\mu$: A function defining the direction of each edge
$\mu:(v_{i},v_{j})\rightarrow[0,1]$
such that for each edge $(v_{i},v_{j})\in E$, $\mu(v_{i},v_{j})=x$ implies $\mu(v_{j},v_{i})=\sqrt{1-x^{2}}$.

In this model, each edge is associated with a scalar $x$ that represents its direction. The value $x$ is a real number in the interval $[0,1]$.
If for an edge $(v_{i},v_{j})$, $\mu(v_{i},v_{j})=x$, then it must hold that $\mu(v_{j},v_{i})=\sqrt{1-x^{2}}$, capturing the edge in both directions.

For example, if $\mu(v_{i},v_{j})=1$ then the edge is an arc (directed edge) connecting node $v_{i}$ to node $v_{j}$. If $\mu(v_{i},v_{j})=0$ then the edge is an arc connecting node $v_{j}$ to node $v_{i}$. If $\mu(v_{i},v_{j})=1/\sqrt{2}$ then the edge is a bidirectional edge connecting node $v_{i}$ to node $v_{j}$ and node $v_{j}$ to node $v_{i}$.

For a scalar-directed edge graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mu)$, the adjacency matrix $\mathbf{A}$ is a square matrix of dimension $|\mathcal{V}|\times|\mathcal{V}|$. The entry $A_{ij}$ of the matrix is defined as follows:

$$ $A_{ij}=\begin{cases}\mu(v_{i},v_{j}),&\text{if }(v_{i},v_{j})\in E\\ 0,&\text{otherwise}\end{cases}$ $$

In this setting, $\mu(v_{i},v_{j})$ captures the direction of the edge from $v_{i}$ to $v_{j}$. It follows that if an edge is present between nodes $v_{i}$ and $v_{j}$ then $A_{ji}=\sqrt{1-A_{ij}^{2}}$.

Therefore, the adjacency matrix captures not only the presence of edges but also their direction according to the function $\mu$.

### D.2 Fuzzy Laplacian Matrix

For a scalar-directed edge graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mu)$ with adjacency matrix $\mathbf{A}$, we define its Fuzzy Laplacian matrix $\mathbf{L}_{F}$ as follows:

The diagonal entries of $\mathbf{L}_{F}$ are zero:

$$ $(\mathbf{L}_{F})_{ii}=0$ $$

The off-diagonal elements are

$$ $(\mathbf{L}_{F})_{ij}=\begin{cases}0&\text{if }A_{ij}=A_{ji}=0\\ e^{i\theta_{ij}}&\text{otherwise}\end{cases}$ (6) $$

where $\theta_{ij}$ is selected such that:

$$ $\cos(\theta_{ij})=A_{ij}$ $$

In other words, the real part of $e^{i\theta_{ij}}$ is equal to the corresponding adjacency matrix entry $A_{ij}$. We require that $0\leq\theta_{ij}\leq\pi/2$. It follows that $\theta_{ji}=\pi/2-\theta_{ij}$.

The Fuzzy Laplacian $\mathbf{L}_{F}$ satisfies the property:

$$ $\mathbf{L}_{F}=i\mathbf{L}_{F}^{*}$ $$

To confirm this, note that $e^{-i\theta_{ij}}=\cos(-\theta_{ij})+i\sin(-\theta_{ij})=\cos(\theta_{ij})-i
\sin(\theta_{ij})$. Therefore, $(\mathbf{L}_{F})_{ji}=\sin(\theta_{ij})+i\cos(\theta_{ij})=ie^{-i\theta_{ij}}$, and thus $\mathbf{L}_{F}=i\mathbf{L}_{F}^{*}$.

The fuzzy Laplacian takes the following form,

$$ $\mathbf{L}_{F}=\begin{pmatrix}0&\cdots&\cdots\\ \vdots&\ddots&e^{i\theta_{ij}}\\ \vdots&ie^{-i\theta_{ij}}&\ddots\\ \end{pmatrix}$ (7) $$

Here, $e^{i\theta_{ij}}$ and $ie^{-i\theta_{ij}}$ are sample off-diagonal elements corresponding to the edge $(v_{i},v_{j})$ in the graph.

### D.3 Properties of Fuzzy Laplacian Matrix 𝐋 F subscript 𝐋 𝐹 \mathbf{L}_{F} bold_L start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT

In this section, we will show that the Fuzzy Laplacian matrix $\mathbf{L}_{F}$ has eigenvalues of the form $a+ia$, where $a\in\mathbb{R}$, and orthogonal eigenvectors.

#### D.3.1 Eigenvalues of the Form a + i ⁢ a 𝑎 𝑖 𝑎 a+ia italic_a + italic_i italic_a

$\mathbf{L}_{F}$ has eigenvalues of the form $a+ia$ with $a\in\mathbb{R}$.

Proof:

Let $\lambda$ be an eigenvalue of $\mathbf{L}_{F}$, and let $\mathbf{w}$ be the corresponding eigenvector. Then:

$$ $\mathbf{L}_{F}\mathbf{w}=\lambda\mathbf{w}\ \ \Rightarrow\ \ \mathbf{w}^{*} \mathbf{L}_{F}^{*}=\lambda^{*}\mathbf{w}^{*}\ \ \Rightarrow\ \ \mathbf{w}^{*} \mathbf{L}_{F}^{*}\mathbf{w}=\lambda^{*}\mathbf{w}^{*}\mathbf{w}$ $$

$$ $\Rightarrow\ \ -i\mathbf{w}^{*}\mathbf{L}_{F}\mathbf{w}=\lambda^{*}\mathbf{w}^ {*}\mathbf{w}\ \ \Rightarrow\ \ -i\lambda\mathbf{w}^{*}\mathbf{w}=\lambda^{*} \mathbf{w}^{*}\mathbf{w}\ \ \Rightarrow\ \ -i\lambda=\lambda^{*}$ $$

where we used $\mathbf{L}_{F}^{*}=-i\mathbf{L}_{F}$ to go to the second line. The last identity holds only when $\lambda=a+ia$ where $a$ is a real number.

#### D.3.2 Orthogonal Eigenvectors

To prove that $\mathbf{L}_{F}$ has orthogonal eigenvectors, we need to show that if $\mathbf{w}$ and $\mathbf{v}$ are eigenvectors corresponding to distinct eigenvalues $\lambda_{1}$ and $\lambda_{2}$ respectively, then $\mathbf{w}$ and $\mathbf{v}$ are orthogonal.

Proof:

Let $\mathbf{L}_{F}\mathbf{w}=\lambda_{1}\mathbf{w}$ and $\mathbf{L}_{F}\mathbf{v}=\lambda_{2}\mathbf{v}$.

$$ $\mathbf{v}^{*}\mathbf{L}_{F}\mathbf{w}=\lambda_{1}\mathbf{v}^{*}\mathbf{w}\ \ \Rightarrow\ \ i\mathbf{v}^{*}\mathbf{L}_{F}^{*}\mathbf{w}=\lambda_{1}\mathbf{ v}^{*}\mathbf{w}\ \ \Rightarrow\ \ i\lambda_{2}^{*}\mathbf{v}^{*}\mathbf{w}= \lambda_{1}\mathbf{v}^{*}\mathbf{w}\ \ \Rightarrow\ \ \lambda_{2}\mathbf{v}^{* }\mathbf{w}=\lambda_{1}\mathbf{v}^{*}\mathbf{w}$ $$

where we used $\lambda=i\lambda^{*}$ derived above for the last step. Therefore, $\lambda_{2}\mathbf{v}^{*}\mathbf{w}=\lambda_{1}\mathbf{v}^{*}\mathbf{w}$.
Since $\lambda_{1}\neq\lambda_{2}$, it must be that $\mathbf{v}^{*}\mathbf{w}=0$, i.e., $\mathbf{w}$ and $\mathbf{v}$ are orthogonal.

Because the eigenvectors of the Fuzzy Laplacian matrix $\mathbf{L}_{F}$ are orthogonal, they can be used as the positional encoding of the nodes of the graph.

## Appendix E Expressivity of Neural Networks using the Fuzzy Laplacian

### E.1 Graph isomorphism for directed graphs with fuzzy edges

First, we extend the standard definition of graph isomorphism to the case of directed graphs with fuzzy edges following a similar approach as in Piperno et al. ().

We only consider graph with a finite number of nodes. Therefore, the set of all edge weights form a countable set with finite cardinality.

Two directed fuzzy graphs $\mathcal{G}=(\mathcal{V}_{\mathcal{G}},\mathcal{E}_{\mathcal{G}},\mu_{\mathcal
{G}})$ and $\mathcal{H}=(\mathcal{V}_{\mathcal{H}},\mathcal{E}_{\mathcal{H}},\mu_{\mathcal
{H}})$ are said to be isomorphic if there exists a bijection $f:\mathcal{V}_{\mathcal{G}}\rightarrow\mathcal{V}_{\mathcal{H}}$ such that, for every pair of vertices $u,v$ in $\mathcal{V}_{\mathcal{G}}$, the following conditions hold:

1. $(u,v)$ is an edge in $\mathcal{E}_{\mathcal{G}}$ if and only if $(f(u),f(v))$ is an edge in $\mathcal{E}_{\mathcal{H}}$ for all vertices $u$ and $v$ in $\mathcal{V}_{\mathcal{G}}$.

2. For every edge $(u,v)$ in $\mathcal{E}_{\mathcal{G}}$ that maps to edge $(f(u),f(v))$ in $\mathcal{E}_{\mathcal{H}}$, the corresponding weights satisfy:

$$ $\mu_{\mathcal{G}}(u,v)=\mu_{\mathcal{H}}(f(u),f(v))$ $$

### E.2 Weisfeiler-Leman test for isomorphism of directed graphs with fuzzy edges

Next, we extend the Weisfeiler-Leman (WL) graph isomorphism test to determine whether two directed graphs with fuzzy edges are isomorphic or not according to the extended definition of graph isomorphism stated above.

WL test is a vertex refinement algorithm that assigns starting features to each node of the graph. The algorithm then aggregates all the features of each node’s neighbors and hashes the aggregated labels alongside the node’s own label into a unique new label. At each iteration of the algorithm, the list of labels are compared across the two graphs. If they are different, the two graphs are not isomorphic. If the labels are no longer updated at each iteration, the two graphs can potentially be isomorphic.

We extend the WL test to directed graphs with fuzzy edges in two ways. Importantly, for all the proofs that follow, we will assume that the graph has a finite number of nodes.

The strong form of WL test for fuzzy directed graphs:

Given a directed graph with fuzzy edges $\mathcal{G}=(\mathcal{V},\mathcal{E},\mu)$, the strong form of the WL test calculates a node coloring $C^{(t)}:\mathcal{V}\rightarrow\{1,2,\ldots,k\}$, a surjective function that maps each vertex to a color. At the first iteration, $C^{(0)}=0$. At all subsequent iterations,

$$ $C^{(t)}(i)=Relabel(C^{{(t-1)}}(i),\{\!\!\{(\mu_{ij},C^{(t-1)}_{j}):j\in N(i)\} \!\!\}),$ (8) $$

where $(\mu_{ij},C^{(t-1)}_{j})$ is tuple comprised of the weight of the edge from node $i$ to node $j$ and $C^{(t-1)}_{j}$, the color of node $j$ in the previous iteration. To simplify notation, we will write $\mu_{ij}$ to denote the weight of the edge connecting nodes $i$ and $j$ instead of $\mu(v_{i},v_{j})$ from now on. $N(i)$ denotes the set of all neighbors of node $i$. Function $Relabel$ is an injective function that assigns a unique new color to each node based on the node’s color in the previous iteration and the tuples formed from its neighbors colors and the value of the edges connecting the node to those neighbors.

From the definition above, it follows that a graph neural network $\Gamma:\mathcal{G}\rightarrow\mathbb{R}^{d}$ that aggregates and updates the node features as follows

$$ $h^{(k)}_{i}=\phi(h^{(k-1)}_{i},f(\{\!\!\{(\mu_{ij},h^{(k-1)}_{j}):j\in N(i)\} \!\!\})),$ (9) $$

will map two graphs $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ to different vectors in $\mathbb{R}^{d}$ if the above strong form of the graph coloring test deems that the two graphs are not isomorphic. In the above equation, $h^{k}_{v}$ is the hidden representation (or feature) of node $i$ at the $k$th layer. $\phi$ and $f$ are injective functions with $f$ acting on multisets of tuples of the edge weights and hidden representation of the neighbors of node $i$.

The proof of above statement is a trivial extension of the proof of Theorem 3 in (Xu et al., ). Briefly, the multiset of the features of the neighbors of node $i$ can be converted from a multiset of tuples of the form $\{\!\!\{(\mu_{ij},h^{(k-1)}_{j}):j\in N(i)\}\!\!\}$ to a multiset of augmented features of the neighbors $\{\!\!\{\tilde{h}^{(k-1)}_{j}:j\in N(i)\}\!\!\}$ where $\tilde{h}^{(k-1)}_{j}=(\mu_{ij},h^{(k-1)}_{j})$. Because $\mu_{ij}$ form a set of finite cardinality, the augmented feature space $\tilde{h}^{(k-1)}_{j}$ is still a countable set and the same proof as in (Xu et al., ) can be applied.

The weak form of WL test for fuzzy directed graphs:

We introduce a weaker form of the WL test for directed graphs with fuzzy edges where a given node cannot distinguish between its neighboring nodes that have the same color based on the weights of the edges that they share. Rather, the node aggregates the weights of all outgoing and incoming edges of its neighbors of a given color.

The weak form of the algorithm is as follows. At the first iteration, $C^{(0)}=0$. At all subsequent iterations,

$$ $\displaystyle C^{(t)}(i)=Relabel$ $\displaystyle\Big{(}C^{(t-1)}(i),$ $\displaystyle\{\big{(}\sum_{j\in N(i)}\delta(C^{(t-1)}_{j}-c)\mu_{ij},$ $\displaystyle\sum_{j\in N(i)}\delta(C^{(t-1)}_{j}-c)\mu_{ji},\ c\big{)}:c\in C ^{(t-1)}_{N(i)}\}\Big{)},$ (10) $$

The term $\sum_{j\in N(i)}\delta(C^{(t-1)}_{j}-c)\mu_{ij}$ is summing over the edge weights $\mu_{ij}$ of all neighbors of $i$ that have color $c$. The tuple now also contains a similar that sums over the edge weights pointing from node $j$ to node $i$, $\mu_{ji}$ for all neighbors of the same color. As stated above, the edge weights are related, namely, $\mu_{ji}=\sqrt{1-\mu_{ij}^{2}}$. The color of the node $i$ at the previous iteration alongside the set of tuples containing the sum of incoming edge weights from neighbors of a given color, the sum of outgoing edge weights to the neighbors of a given color, and the color of those neighbors are inputs to the injective function $Relabel$, which assigns a unique color to node $i$ for the next iteration.

Theorem. Let $\Gamma:\mathcal{G}\rightarrow\mathbb{R}^{d}$ be a graph neural network that updates node features as follows,

$$ $h^{(k)}_{i}=MLP^{(k)}\left(h^{(k-1)}_{i},\Re\left[\sum_{j\in N(i)}(\mathbf{L}_ {F})_{ij}h^{(k-1)}_{j}\right],\Im\left[\sum_{j\in N(i)}(\mathbf{L}_{F})_{ij}h^ {(k-1)}_{j}\right]\right),$ (11) $$

where $MLP^{(k)}$ is a multi-layer perceptron. The last two terms of the MLP input are the real and imaginary parts of the aggregated features of the neighbors of node $i$, $\sum_{j\in N(i)}(\mathbf{L}_{F})_{ij}h^{(k-1)}_{j}$. With sufficient number of layers, the parameters of $\Gamma$ can be learned such that it is as expressive as the weak form of the WL test, in that $\Gamma$ maps two graphs $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ that the weak form of WL test decides to be non-isomorphic to different embeddings in $\mathbb{R}^{d}$.

Proof.
First, we claim that there exists a GNN of the form:

$$ $h^{(k)}_{i}=\phi\left(f(h^{(k-1)}_{i}),\sum_{j\in N(i)}\mu_{ij}f(h^{(k-1)}_{j} ),\sum_{j\in N(i)}\mu_{ji}f(h^{(k-1)}_{j})\right),$ (12) $$

that is as expressive as the weak form of the WL test. We will prove this by induction.

Note that $F(i)=\sum_{j\in N(i)}\mu_{ij}f(h^{(k-1)}_{j})$ and $G(i)=\sum_{j\in N(i)}\mu_{ji}f(h^{(k-1)}_{j})$ are only injective up to the sum of the incoming and outgoing weights for a given type of neighbor (see Figure E.2 for an example). This is because the weak form of the WL test only accounts for the sum of the weights of outgoing and incoming edges to neighboring nodes of a given color. This is different from previous work that dealt with undirected graphs (Xu et al., ) or graphs with directed but unweighted edges (Rossi et al., ).

Figure: Figure E.2: Examples of two neighborhoods of a node $i$ are shown on the left and right. The strong form of the WL test for directed graphs with fuzzy edges can distinguish these two neighborhoods. The weak form of the WL test, however, cannot do so because the sum of the weights of the edges connecting node $i$ to green-colored neighbors is 0.9 in both cases.
Refer to caption: extracted/5934640/figs/figA2.png

Let’s denote with $X$ the multiset of the colors of all neighbors of node $i$ after a given number of iterations of the color refinement algorithm. Equivalently, we can consider the multiset of all the features of the neighboring nodes of $i$ after a given number of iteration of the GNN algorithm.

Because our graphs have a finite number of nodes, $X$ is bounded. For graphs without directed edges or weighted edges, it can be shown (see Lemma 5 of (Xu et al., )) that there always exists an injective function $f$ such that $F(x)=\sum_{x\in X}f(x)$ is injective for all multisets $X$, i.e. if $F(X)=F(Y)$ for two multisets $X$ and $Y$ then $X=Y$. This result can be easily extended to regular directed graphs without weighted edges by considering the neighbors with incoming edges separately from the neighbors with outgoing edges (Rossi et al., ). For such graphs, a function $f$ exists such that the tuple $(\sum_{x\in X_{\rightarrow}}f(x),\sum_{x\in X_{\leftarrow}}f(x))$ is only the same for two nodes if they have identical neighborhoods. $X_{\rightarrow}$ and $X_{\leftarrow}$ denotes the multiset of features of the neighboring nodes that are connected to $i$ using an outgoing edge and incoming edge respectively.

Let’s proceed with our proof by induction. At the first iteration, $k=0$, all the nodes have the same color corresponding to the same trivial feature (e.g. scalar 0). Functions $F(i)$ and $G(i)$ then simply sum the weights of the outgoing and incoming edges of node $i$ respectively and use their values to assign updated feature $h_{i}^{(1)}$. This is identical to the procedure that the weak WL algorithm is using to assign new colors to the nodes at its first iteration. Therefore, nodes that would be assigned a given color under the WL color refinement algorithm at its first iteration will also be assigned the same updated feature vector $h_{i}^{(1)}$.

Assume that our claim hold for iteration $k-1$. This mean that all the nodes that the weak WL color refinement assigns to different colors are also assigned different features $h_{i}^{(k-1)}$ by the GNN. Following the weak WL test, at iteration $k$, we need to sum the incoming and outgoing edges across all the neighbors with the same $h_{j}^{(k-1)}$. An example of an $f$ that allows this is one-hot encoding of all features. There always exists a number $N\in\mathbb{N}$ such that the features $h_{i}^{(k-1)}$ across all nodes $i$ of the graph can be one-hot encoded in an $N$ dimensional vector. It follows,

$$ $\sum_{j\in N(i)}\mu_{ij}f(h^{(k-1)}_{j})=\left[\sum_{j\in N(i)}\delta(h^{(k-1) }_{j}-h)\mu_{ij}\right]_{h\in H^{(k-1)}_{N(i)}},$ (13) $$

where $H^{(k-1)}_{N(i)}$ is the set of all the features $h_{j}^{(k-1)}$ of nodes $j$ that are neighbors of node $i$. A similar expression can be written for the sum over the incoming weights $\mu_{ji}$. Thus, at iteration $k$ the GNN will assign distinct features to all nodes that are also assigned a distinct color at iteration $k$ of the weak WL color refinement algorithm.

To complete the proof, we note that because $\mu_{ji}=\sqrt{1-\mu_{ij}^{2}}$, we can reparameterize the edge weights as $\mu_{ij}=\cos(\theta_{ij})$ and $\mu_{ji}=\sin(\theta_{ij})$ with $0\leq\theta_{ij}\leq\pi/2$. Eq. 12 can be rewritten as:

$$ $h^{(k)}_{i}=\phi\left(f(h^{(k-1)}_{i}),\Re\left[\sum_{j\in N(i)}(\mathbf{L}_{F })_{ij}f(h^{(k-1)}_{j})\right],\Im\left[\sum_{j\in N(i)}(\mathbf{L}_{F})_{ij}f (h^{(k-1)}_{j})\right]\right),$ (14) $$

Finally, we use universal approximation theorem of multilayer perceptrons (Hornik et al., ) to model both the $\phi$ computation at layer $k$ and the $f$ computation for the next layer, $k+1$. Namely $MLP^{(k)}$ denotes $f^{(k+1)}\circ\phi^{(k)}$. Taken together, the GNN in Eq. 11 is as expressive as the weak form of the WL color refinement algorithm for directed graphs with fuzzy edges.

## Appendix F Magnetic Laplacian is not as expressive as the fuzzy Laplacian.

### F.1 Magnetic Laplacian

A commonly used Laplacian for directed graphs is the Magnetic Laplacian. To construct the Magnetic Laplacian, we start with the asymmetric adjacency matrix of the graph and symmetrize it. Note that we are using the conventional definition of adjacency matrix in this section and not the fuzzy definition used above. The adjacency matrix $\mathbf{A}$ for a directed graph is a binary matrix, where $A_{ij}=1$ represents the presence of a directed edge from vertex $v_{i}$ to vertex $v_{j}$ and $A_{ij}=0$ represents the absence of such an edge.

The symmetrized adjacency matrix $S$ is defined as:

$$ $S_{ij}=\frac{1}{2}(A_{ij}+A_{ji}).$ $$

To capture the direction of the edges, a phase matrix is defined as,

$$ $\Theta^{(q)}_{ij}=2\pi q(A_{ij}-A_{ji}).$ (15) $$

The Hermitian adjacency matrix is defined as the element-wise product of the above two matrices,

$$ $\mathbf{H}^{(q)}=S\odot\Theta({q}).$ $$

$\mathbf{H}^{(q)}$ has some useful properties in capturing the directionality of the edges of the graph. For example, for $q=1/4$, if there is an edge connecting $j$ to $k$ but no edge connecting $k$ to $j$ then $\mathbf{H}^{(q)}_{jk}=i/2$ and $\mathbf{H}^{(q)}_{kj}=-i/2$. Although this encoding of edge direction is useful, the fact that both incoming and outgoing edges are purely imaginary (only different up to a sign) means that in general it is impossible to distinguish features aggregated from neighbors that are connected to a node through outgoing edges from those connected through incoming edges. The fuzzy Laplacian $\textbf{L}_{F}$, however, trivially distinguishes features from outgoing neighbors from those of incoming neighbors by keeping one set real and the other imaginary. We will expand on the implication of this for the expressivity of graph neural networks constructed using these two approaches below.

### F.2 Limitations of the magnetic Laplacian

Conventionally, a magnetic Laplacian is defined by subtracting $\mathbf{H}^{(q)}$ from the degree matrix, $\mathbf{L}_{M}=\mathbf{D}-\mathbf{H}^{(q)}$. If such a Laplacian matrix is used to aggregate information of the nodes of the graph, the features of a node itself are combined with the features of its neighboring nodes. For directed graphs, there are two categories of neighboring nodes, those connected to a node $i$ with outgoing edges from $i$ ($A_{ij}=1$ and $A_{ji}=0$) and those connected to node $i$ with incoming edges ($A_{ji}=1$ and $A_{ij}=0$). Of course, it is possible for a neighboring node to be both an outgoing and incoming neighbor, in which case $A_{ji}=A_{ij}=1$.

Features from these two categories of neighbors, alongside the feature of the node itself, form three distinct categories that in general must be kept distinct for maximum expressivity (Rossi et al., ). Using a complex Laplacian for aggregating the features across the nodes allows a simple mechanism for distinguishing two categories (corresponding to the real and imaginary parts of the resulting complex features) but not all three. To get around this limitation, we keep the self-feature of every node distinct from the aggregated features of its neighbors and concatenate it with the aggregated feature prior to applying the multi-layer perceptron to update the features from one layer to the next. Therefore, to maximize the expressivity of the magnetic Laplacian, we set the diagonal terms of the Laplacian matrix to zero, and define the magnetic Laplacian simply as $\mathbf{L}_{M}=\mathbf{H}^{(q)}$.

Lemma. For simple directed graphs without fuzzy edges, using the magnetic Laplacian $\mathbf{L}_{M}$ in the graph neural network of Eq. 11 in the place of the fuzzy Laplacian $\textbf{L}_{F}$,

$$ $h^{(k)}_{i}=MLP^{(k)}\left(h^{(k-1)}_{i},\Re\left[\sum_{j\in N(i)}(\mathbf{L}_ {M})_{ij}h^{(k-1)}_{j}\right],\Im\left[\sum_{j\in N(i)}(\mathbf{L}_{M})_{ij}h^ {(k-1)}_{j}\right]\right),$ (16) $$

does not decrease the expressivity of the graph neural network in that both networks are as expressive as the weak form of the WL graph isomorphism test.

Proof.. To prove this statement, we will show that the MLP can learn a linear combination of its input that would it make it equivalent to the graph neural network defined in Eq. 11. This is because when $q=1/8$, the real and imaginary parts of the aggregated features of the neighboring nodes are linear combinations of the features of the outgoing and incoming neighbors.

The fuzzy Laplacian $\textbf{L}_{F}$ directly separates the features of outgoing and incoming neighbors in the real and imaginary components of the aggregated features respectively. Define $F_{\rightarrow}(i)=\Re\left[\sum_{j\in N(i)}(\mathbf{L}_{F})_{ij}h^{(k-1)}_{j}\right]$ and $F_{\leftarrow}(i)=\Im\left[\sum_{j\in N(i)}(\mathbf{L}_{F})_{ij}h^{(k-1)}_{j}\right]$ as the features aggregated from the outgoing and incoming neighbors respectively. With $q=1/8$, we have,

$$ $\displaystyle\Re\left[\sum_{j}(\mathbf{L}_{M})_{ij}h^{(k-1)}_{j}\right]$ $\displaystyle=\frac{1}{2\sqrt{2}}(F_{\rightarrow}(i)+F_{\leftarrow}(i))$ $\displaystyle\Im\left[\sum_{j}(\mathbf{L}_{M})_{ij}h^{(k-1)}_{j}\right]$ $\displaystyle=\frac{1}{2\sqrt{2}}(F_{\rightarrow}(i)-F_{\leftarrow}(i))$ $$

The MLP layer then in the graph neural network of Eq.16 simply needs to learn a linear combination of its second and third concatenated inputs to become equivalent to the graph neural network of Eq.11. Therefore, the two graph neural networks are equally as expressive.

Next, we extend the definition of the magnetic Laplacian to directed graphs with fuzzy edges, defined above, which have weighted edges that indicate intermediate values of directionality between the two extremes of the edge pointing form node $i$ to node $j$ and the edge pointing from node $j$ to node $i$. We implemented these intermediate values in the fuzzy Laplacian by assigning a weight $\mu_{ij}$ to each edge. Although not necessary in general, we further assumed that $\mu_{ji}=\sqrt{1-\mu_{ij}^{2}}$. This constraint allowed us to think of each edge weight as an angle $\theta_{ij}=\arccos(\mu_{ij})$ which we then used to define the fuzzy Laplacian matrix (Eq.7).

To extend the magnetic Laplacian to directed graphs with fuzzy edges, we can assign a different value for $q$ (Eq. 15) to each edge. To keep the notation comparable to that of the fuzzy Laplacian, we instead assign a separate angle $\theta_{ij}$ to each edge and define the magnetic Laplacian as, $(\mathbf{L}_{M})_{ij}=e^{i\theta_{ij}}$ for $j>i$ if there is an edge between the nodes $i$ and $j$, $S_{ij}\neq 0$. To ensure that magnetic Laplacian remains Hermitian, $\mathbf{L}_{M}^{*}=\mathbf{L}_{M}$, we require that $(\mathbf{L}_{M})_{ji}=e^{-i\theta_{ij}}$, which defines the $\theta_{ij}$ values for $j<i$. It follows that $\theta_{ji}=-\theta_{ij}$.

$$ $\mathbf{L}_{M}=\begin{pmatrix}0&\cdots&\cdots\\ \vdots&\ddots&e^{i\theta_{ij}}\\ \vdots&e^{-i\theta_{ij}}&\ddots\\ \end{pmatrix}$ $$

It remains to be shown how we can assign the $\theta_{ij}$ values from the $\mu_{ij}$ in a self-consistent way.

Lemma. For directed graphs with fuzzy edges, the graph neural network constructed using the magnetic Laplacian $\mathbf{L}_{M}$, Eq. 16, is not as expressive as the weak form of the WL graph isomorphism test, and therefore not as expressive as the graph neural network defined using the fuzzy Laplacian $\mathbf{L}_{F}$, Eq. 11.

Proof. The Fuzzy Laplacian matrix conveniently captures the outgoing and incoming weights of the edges from one node to another, $(\mathbf{L}_{F})_{ij}=\cos(\theta_{ij})+i\sin(\theta_{ij})$. Namely, $\Re(\mathbf{L}_{F})_{ij}=\cos(\theta_{ij})$ is the outgoing weight $\mu_{ij}$ of node $i$ to node $j$ and $\Im(\mathbf{L}_{F})_{ij}=\sin(\theta_{ij})$ is in the incoming weight from node $j$ to node $i$, $\mu_{ji}$. Importantly, the outgoing weight from node $i$ to node $j$, $\mu_{ij}$, is the same as the incoming weight from node $i$ to node $j$. Similarly, the incoming weight from node $j$ to node $i$, $\mu_{ji}$ is the same as the outgoing weight from node $j$ to node $i$. These relationships are directly captured in the fuzzy Laplacian because $(\mathbf{L}_{F})_{ji}=\sin(\theta_{ij})+i\cos(\theta_{ij})$.

Let’s consider how we can construct the magnetic Laplacian $\mathbf{L}_{M}$ from $\mu_{ij}$. The weight of the edge from $i$ to $j$ is $\mu_{ij}$. The weight from $j$ to $i$ is $\mu_{ij}=\sqrt{1-\mu_{ij}^{2}}$.

From our definition above, $(\mathbf{L}_{M})_{ij}=\cos(\theta_{ij})+i\sin(\theta_{ij})$ and $(\mathbf{L}_{M})_{ji}=\cos(\theta_{ij})-i\sin(\theta_{ij})$. To relate $\theta_{ij}$ to $\mu_{ij}$ we note that the ratio of outgoing and incoming weights for node $i$ is $\mu_{ij}/\mu_{ji}$ and for node $j$ is $\mu_{ji}/\mu_{ij}$. Because these two quantities are reciprocals of each other, we can define

$$ $\ln\frac{\mu_{ij}}{\mu_{ji}}=\tan(2\theta_{ij}).$ $$

We chose the $\tan$ function on the right hand side because $\tan(-2\theta_{ij})=-\tan(2\theta_{ij})$. Moreover, the log-ratio on the left hand side of above equation can take on any value from $-\infty$ to $\infty$. The right hand side can take on a similar range of values if we allow $-\pi/4\leq\theta_{ij}\leq\pi/4$. The specific form of this function does not matter as long it is an odd function that spans the specified domain and range.

Assume that each node of the graph has a trivial scalar feature equal to 1 (the first iteration of the WL test). The magnetic Laplacian applied to this graph to aggregate the features of the neighbors of node $i$ gives,

$$ $\displaystyle\Re\left[\sum_{j}(\mathbf{L}_{M})_{ij}\right]$ $\displaystyle=\sum_{j}\cos\left(\frac{1}{2}\arctan\left[\ln\frac{\mu_{ij}}{ \sqrt{1-\mu_{ij}^{2}}}\right]\right)$ $\displaystyle\Im\left[\sum_{j}(\mathbf{L}_{M})_{ij}\right]$ $\displaystyle=\sum_{j}\sin\left(\frac{1}{2}\arctan\left[\ln\frac{\mu_{ij}}{ \sqrt{1-\mu_{ij}^{2}}}\right]\right)$ $$

The above functions for a single value of $\mu_{ij}$ are plotted in Figure F.3. Importantly, these functions are not an injective function of all possible neighborhoods of node $i$. Consider the case of a node $i$ that has 5 neighbors with $\mu_{ij}$ values that result in $\cos(\theta_{ij})=0.8$ and $\sin(\theta_{ij})=0.6$, and another set of 5 neighbors with $\mu_{ij}$ values that result in $\cos(\theta_{ij})=0.8$ and $\sin(\theta_{ij})=-0.6$. It follows that $\Re\left[\sum_{j}(\mathbf{L}_{M})_{ij}\right]=8$ and $\Im\left[\sum_{j}(\mathbf{L}_{M})_{ij}\right]=0$. The same values would have been obtained if node $i$ had a neighborhood comprised of 8 neighbors with $\mu_{ij}$ values such that $\cos(\theta_{ij})=1$ and $\sin(\theta_{ij})=0$. Therefore, the magnetic Laplacian cannot distinguish distinct neighborhoods and is not as expressive as the weak form of the WL test.

Figure: Figure F.3: The functions used to map the edge weights $\mu_{ij}$ to the $\theta_{ij}$ values of the magnetic Laplacian. $\cos(\theta_{ij})$ is plotted on the left. $\sin(\theta_{ij})$ is plotted in the middle. $\sin(\theta_{ij})$ versus $\cos(\theta_{ij})$ is plotted on the right.
Refer to caption: extracted/5934640/figs/figA3.png

In practice, the more important limitation of the magnetic Laplacian is that it aggregates neighborhood information that results in linear combinations of the outgoing and incoming features to a node. In contrast, the fuzzy Laplacian by construct always keeps these contributions separate in the real and imaginary parts of the aggregated feature vector. In many applications, we would like to access the self-feature of the node and the incoming and outgoing aggregated features of its neighbors separately. The fuzzy Laplacian is thus a better choice. It could be argued that the MLP of a graph neural network can learn to disentangle the linear combinations of the outgoing and incoming features aggregated by the magnetic Laplacian. In general, however, this is not possible by the types of the graph neural networks that we considered here. This is because the parameters are MLP are the same for all the nodes of the graph. The linear combinations, however, depend on the specific weights $\mu_{ij}$ connecting node $i$ to its neighboring nodes $j$. Therefore, in general, the MLP will not be able to learn to disentangle the linear combinations of the incoming and outgoing features aggregated by the magnetic Laplacian.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: 1. What is the runtime of the proposed methods, and how does it compare to the baselines? Providing a detailed runtime analysis would enhance the evaluation.
2. For the node classification task, how does the performance change if the edge direction is made learnable? Prior work, such as [3], suggests that it will likely improve performance.
soundness: 3
strengths: The paper is well-written and offers an intriguing evaluation on graph ensemble data, where the graph structure remains fixed but multiple realizations of node features and targets exist—an important task that demands attention. Additionally, the authors demonstrate that their method is more expressive than the Magnet Laplacian.
summary: Traditional GNNs often suffer from feature oversmoothing due to their diffusion-like message-passing mechanism over undirected graphs. To address this, the authors propose the Continuous Edge Direction (CoED) GNN framework that assigns differentiable fuzzy edge directions to graph edges, allowing features to preferentially flow in one direction and introducing a complex-valued Laplacian for directed graphs with fuzzy edges.
weaknesses: 1. The paper omits important baselines and related work, such as [3], which share similar concepts and should be included for a comprehensive comparison.
2. If the theorem is highlighted in the abstract and introduction, it should be mentioned explicitly in the main paper for consistency and emphasis.
3. The evaluation on node classification tasks is unconvincing:
   - Only one homophilous dataset is included in the study.
   - The Chameleon and Squirrel datasets used are known to be problematic [1], and it is advisable to use filtered versions for more reliable results.
   - The heterophilous datasets used for evaluation are all quite small, and the results do not approach the current state-of-the-art. To more robustly validate the claims, experiments on larger-scale heterophilous datasets [1][2] are necessary.


References

[1] Platonov, O., Kuznedelev, D., Diskin, M., Babenko, A., & Prokhorenkova, L. (2023). A critical look at the evaluation of GNNs under heterophily: Are we really making progress?. arXiv preprint arXiv:2302.11640.

[2] Derek Lim, Felix Hohne, Xiuyu Li, Sijia Linda Huang, Vaishnavi Gupta, Omkar Bhalerao, and Ser Nam Lim. Large scale learning on non-homophilous graphs: New benchmarks and strong simple methods. Advances in Neural Information Processing Systems, 34, 2021.

[3]Finkelshtein, Ben, et al. "Cooperative graph neural networks." ICML 2024.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 3
questions: See Weaknesses
soundness: 3
strengths: 1. The paper is well-written, and the motivation of the model is reasonable.
2. The idea of learning optimal edge directions, and using phase angles to scale the direction-aware information flow is intuitive and interesting.
3. Theoretical Analysis appears to be rigorous and formally describes the expressive power of the proposed model
summary: The paper presents a novel graph neural network architecture named Continuous Edge Direction (CoED) GNN, which improves upon traditional GNNs by introducing the concept of fuzzy edge directions. These directions allow for more effective feature propagation by varying continuously, enabling more controlled information flow between nodes. CoED GNN utilizes a complex-valued Laplacian that represents these fuzzy directions, where the real and imaginary parts correspond to opposite directions of information flow. The architecture is especially effective in graph ensemble settings, where the structure is static but node features vary, such as in gene regulatory networks and power grids.
weaknesses: 1. Directed graph Laplacians have been extensively explored in previous studies, including seminal works like those by F. Chung [1] and on Hermitian Laplacian [2]. The introduction of the fuzzy graph Laplacian is a significant contribution of this work. However, while comparisons have been made with the magnetic Laplacian, it would be beneficial to extend these comparisons to include contemporary directed graph Laplacians to fully contextualize its advancements and differences.

2. From my perspective, employing phase angles to dictate continuous edge directions is designed to adaptively modulate the weights from in-neighbors and out-neighbors. However, could we alternatively employ a shared learnable function $f$, , such as, $\mathrm{Sigmoid}(f(u_i, u_j))$ and $1- \mathrm{Sigmoid}(f(u_i, u_j))$ , to scale the information flow? This approach might offer a simpler alternative to using complex values. Furthermore, related work such as superGAT [3] enhances GAT by identifying important neighbors, a concept that could potentially be extended to directed graphs. By incorporating the separate message-passing strategy of DirGNN—which independently considers in-edges and out-edges—we could apply superGAT to each direction separately. Subsequently, the two representations could be combined to produce the final node representation.

3. While the concept of learning optimal edge directions is intriguing, I find the experimental results presented in the paper unconvincing. Firstly, as illustrated in Table 1, the CoED model achieves the best results in only one out of the five datasets examined. Secondly, the paper omits experiments on several pertinent datasets commonly used in related works, notably those employed by DirGNN, which is frequently referenced in this study. Specifically, there are no results included for datasets such as SNAP-PATENTS, ROMAN-EMPIRE, and ARXIV-YEAR, despite DirGNN reporting significant findings on these along with various other baselines. Incorporating these datasets could enhance the persuasiveness of the results. Thirdly, I view SAGE as a specific instantiation of DirGNN, implying that the performance of DirGNN should be at least as good, if not better, than that of SAGE. To enhance DirGNN's performance, particularly on the Texas and Wisconsin datasets, I experimented with incorporating an MLP layer following the final layer of DirGNN. This modification has significantly improved its performance. Additionally, the Texas and Wisconsin datasets are relatively small, and the Chameleon and Squirrel datasets reportedly contain significant amounts of duplicate data [4]. Therefore, I believe these datasets may not be the most suitable for evaluating the model's performance. The authors might consider using the AM-Photo, AM-Computer, or Citeseer datasets, as employed in [5], to provide a more robust assessment.

   

[1] Laplacians and the Cheeger inequality for directed graphs. 

[2] Graph signal processing for directed graphs based on the hermitian laplacian.

[3] How to find your friendly neighborhood: Graph attention design with self-supervision. 

[4] A critical look at the evaluation of GNNs under heterophily: Are we really making progress?

[5] Digraph Inception Convolutional Networks.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. Could the authors provide theoretical analysis of the superiority of CoED GNN over GAT?
2. In Equation (1), why $(L_F)_{ij}=\text{exp}(i\theta_{ij})$? Is the first $i$ the imaginary unit? If so, the authors should use different letter to denote the $i$th node.
3. Could the authors provide explanation and experimental results to show that CoED GNN can alleviate the problem of homogenization of information across the graph?
soundness: 3
strengths: 1. The fuzzy Laplacian enables CoED GNN to handle directional message passing, preserving the discriminative features over longer distances.
2. CoED GNN shows performance gains in diverse datasets, highlighting its adaptability to both synthetic and real-world applications.
summary: This paper leverages "fuzzy edges" to improve the flow of information across a graph. CoED GNN learns directions (coefficients) for each edge. This is achieved through a complex-valued fuzzy Laplacian. Experiments on both synthetic and real-world datasets demonstrate CoED GNN's good performance in node classification and regression tasks, attributed to its ability to optimize long-range feature propagation.
weaknesses: 1. The edge directions are essentially edge weights. GAT can also learn edge weights, and it even leverages node features while doing so, unlike CoED GNN. Since GAT utilizes more information, it should perform better. At the very least, the weights learned by GAT should not be worse than those learned by CoED GNN. It is puzzling that GAT performs worse than CoED GNN in the experimental results provided by the authors. The authors should provide theoretical provements to show that CoED GNN is superior.
2. The authors should explain why CoED GNN can alleviate the problem of homogenization of information across the graph. Also, the authors should  provide experimental validations.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: - What is the runtime comparison explicitly, e.g. via a table, when comparing CoDE with MagNet etc.?
- The paper uses the idea of separate aggregation of in- and out- information in addition to self. [1] does something similar and is based on flow imbalance in its objective. How are they compared empirically?
- In line 315, where should $\alpha$ appear in the formula?
- Table 1 compares a CoED version without edge direction learning, how about with edge direction learning for node classification?

Reference:

[1] He, Y., Reinert, G., & Cucuringu, M. (2022, December). DIGRAC: digraph clustering based on flow imbalance. In Learning on Graphs Conference (pp. 21-1). PMLR.
soundness: 3
strengths: - The paper is generally well-written and well-motivated.
- The topics of expressive power, directed GNNs, and oversmoothing in GNNs are significant.
- The proposal of fuzzy edge directions is novel and effective.
- The statements are theoretically and empirically grounded.
summary: This paper proposes a novel framework for undirected and directed graphs called Continuous Edge Direction (CoED) GNN, based on the insight to assign fuzzy edge directions in a novel complex-valued Laplacian matrix. The paper includes empirical comparison to various GNNs on both synthetic and real-world datasets, and theoretically shows the expressive power of CoED.
weaknesses: - No source code is provided.
- It would be better to show the runtime comparison explicitly, e.g. via a table, when comparing CoDE with MagNet etc.
- The paper uses the idea of separate aggregation of in- and out- information in addition to self. [1] does something similar and is based on flow imbalance in its objective. It would be better to add some discussion and comparison in experiments.
- In line 315, it is not clear where $\alpha$ should appear in the formula.
- Table 1 compares a CoED version without edge direction learning, how about with edge direction learning for node classification?
- (bonus) It would be better to show the ultimate directions (visualized) for some directed or undirected graphs to show more insights.

Reference:

[1] He, Y., Reinert, G., & Cucuringu, M. (2022, December). DIGRAC: digraph clustering based on flow imbalance. In Learning on Graphs Conference (pp. 21-1). PMLR.

Meta-review: The paper presents a novel framework for graph neural networks (GNNs) that introduces fuzzy edge directions using a complex-valued Laplacian. This approach addresses key challenges in traditional GNNs, such as feature oversmoothing and limited expressiveness for directed graphs. The proposed CoED GNN enhances information flow between nodes by assigning differentiable edge directions, enabling controlled and adaptive message passing. This framework shows promise, particularly in scenarios involving static graph structures but varying node features, such as graph ensemble tasks. The authors provide a solid theoretical foundation for the expressiveness of their model and demonstrate its effectiveness on both synthetic and real-world datasets.

**Strengths**

Novelty: The introduction of fuzzy edge directions and a complex-valued Laplacian is innovative and provides a fresh perspective on improving GNN performance.

Empirical Evaluation: The evaluation includes both synthetic and real-world datasets, showcasing the model's adaptability across diverse scenarios.

Relevance to Graph Ensemble Tasks: The paper addresses graph ensemble tasks, a valuable yet underexplored application area in graph learning.

**Weaknesses**
While the reviewers raised concerns regarding the selection of datasets (especially commonly used heterophilic ones) and baseline comparisons, I find these points less critical given the primary application focus on graph ensemble tasks. However, I identified a few areas for improvement that could strengthen the final version:

Potential Overfitting and Task Limitations. Assigning a freely trainable parameter to each edge introduces the risk of overfitting, particularly in tasks outside the graph ensemble domain where the number of parameters might far exceed the available labels. This concern aligns with the authors' response to Reviewer b7qy, where masking nodes resulted in notable performance degradation. I suggest to emphasize the primary focus on graph ensemble tasks in the introduction/motivation and explicitly discuss the model's limitations on commonly used heterophilic graphs. This clarification would help future researchers better understand the scope of the approach.

Expressive Power and WL Test Comparisons. The proposed framework's expressive power comparisons with standard GNNs and magnetic Laplacians, framed through the lens of WL tests, require clarification. The framework's edge-specific parameters inherently break permutation symmetry, which the WL tests assume to preserve. I suggest to address this distinction explicitly in the final version and provide a more extensive comparison with magnetic Laplacians, particularly regarding expressive power. Additionally, several relevant works on magnetic Laplacians were not discussed. I recommend incorporating those fundamental references to provide a more comprehensive discussion:

1. MA Shubin. Discrete Magnetic Laplacian. Communications in Mathematical Physics, 1994.

2. Zhang et al. Magnet: A Neural Network for Directed Graphs. NeurIPS 2021.

3. Fanuel et al. Magnetic Eigenmaps for Community Detection in Directed Networks. Physical Review E, 2017.

4. Fanuel et al. Magnetic Eigenmaps for the Visualization of Directed Networks. Applied and Computational Harmonic Analysis, 2018.

5. Furutani et al. Graph Signal Processing for Directed Graphs Based on the Hermitian Laplacian. ECML PKDD 2019.

6. Geisler et al. Transformers Meet Directed Graphs. ICML 2023.

7. Huang et al. What Are Good Positional Encodings for Directed Graphs? arXiv 2024.


Despite the identified weaknesses, I believe the paper's strengths in novelty and relevance outweigh its limitations. The innovative contributions and targeted applications suggest a promising direction for future research, leading me to lean toward acceptance.

Decision: Accept (Poster)