From Hypergraph Energy Functions to Hypergraph Neural Networks
Hypergraphs are a powerful abstraction for representing higher-order
interactions between entities of interest. To exploit these relationships in
making downstream predictions, a variety of hypergraph neural network
architectures have recently been proposed, in large part building upon
precursors from the more traditional graph neural network (GNN) literature.
Somewhat differently, in this paper we begin by presenting an expressive family
of parameterized, hypergraph-regularized energy functions. We then demonstrate
how minimizers of these energies effectively serve as node embeddings that,
when paired with a parameterized classifier, can be trained end-to-end via a
supervised bilevel optimization process. Later, we draw parallels between the
implicit architecture of the predictive models emerging from the proposed
bilevel hypergraph optimization, and existing GNN architectures in common use.
Empirically, we demonstrate state-of-the-art results on various hypergraph node
classification benchmarks. Code is available at
https://github.com/yxzwang/PhenomNN.
Yuxin Wang, Quan Gan, Xipeng Qiu, Xuanjing Huang, David Wipf
# From Hypergraph Energy Functions to Hypergraph Neural Networks

###### Abstract

Hypergraphs are a powerful abstraction for representing higher-order interactions between entities of interest. To exploit these relationships in making downstream predictions, a variety of hypergraph neural network architectures have recently been proposed, in large part building upon precursors from the more traditional graph neural network (GNN) literature. Somewhat differently, in this paper we begin by presenting an expressive family of parameterized, hypergraph-regularized energy functions. We then demonstrate how minimizers of these energies effectively serve as node embeddings that, when paired with a parameterized classifier, can be trained end-to-end via a supervised bilevel optimization process. Later, we draw parallels between the implicit architecture of the predictive models emerging from the proposed bilevel hypergraph optimization, and existing GNN architectures in common use. Empirically, we demonstrate state-of-the-art results on various hypergraph node classification benchmarks. Code is available at [https://github.com/yxzwang/PhenomNN](https://github.com/yxzwang/PhenomNN).

Machine Learning, Hypergraph Energy Functions, Hypergraph Neural Networks

## 1 Introduction

Hypergraphs represent a natural extension of graphs, whereby each hyperedge can link an arbitrary number of hypernodes (or nodes for short). This flexibility more directly facilitates the modeling of higher-order relationships between entities (Chien et al., 2022; Benson et al., 2016, 2017) leading to strong performance in diverse real-world situations (Agarwal et al., 2005; Li and Milenkovic, 2017; Feng et al., 2019; Huang and Yang, 2021). Currently, hypergraph-graph-based modeling techniques frequently rely, either implicitly or explicitly, on some type of expansion (e.g., clique, star), which effectively converts the hypergraph into a regular graph with a new edge set and possibly additional nodes as well. For example, one approach is to first extract a particular expansion graph and then build a graph neural network (GNN) model on top of it (Zhang et al., 2022).

We instead adopt a different starting point that both allows us to incorporate multiple expansions if needed, but also transparently explore the integrated role of each expansion within a unified framework. To accomplish this, our high-level strategy is to first define a family of parameterized hypergraph energy functions, with regularization factors that we later show closely align with popular existing expansions. We then demonstrate how the minimizers of such energy functions can be treated as learnable node embeddings and trained end-to-end via a bilevel optimization process. Namely, the lower-level minimization process produces optimal features contigent on a given set of parameters, while the higher-level process trains these parameters (and hence the features they influence) w.r.t. downstream node classification tasks.

To actualize this goal, after presenting related work in Section 2, we provide relevant background and notation w.r.t. hypergraphs in Section 3. The remainder of the paper then presents our primary contributions, which can be summarized as follows:

* We present a general class of hypergraph-regularized energy functions in Section 4 and elucidate their relationship with traditional hypergraph expansions that have been previously derived from spectral graph theory.
* We demonstrate how minimizers of these energy functions can serve as principled, trainable features for hypergraph prediction tasks in Sections 5 and 6. And by approximating the energy minimizers using provably-convergence proximal gradient steps, the resulting architecture borrows the same basic structure as certain graph neural network layers that: (i) have been fine-tuned to accommodate hypergraphs, and (ii) maintain
the inductive bias infused by the original energy function.
* The resulting framework, which we name **PhenomNN** for _Purposeful Hyper-Edges iN Optimization Motivated Neural Networks_, is applied to a multitude of hypergraph node classification benchmarks in Section 7, achieving competitive or SOTA performance in each case.

## 2 Related Work

**Hypergraph Expansions/Neural Networks**. Hypergraphs are frequently transformed into graphs by expansion methods including the clique and star expansions. An extensive spectral analysis study of of different hypergraph expansions is provided in (Agarwal et al., 2006), but not from the vantage point of energy functions as is our focus. An alternative line expansion (Yang et al., 2020) has also been proposed that can be viewed in some sense as a hybrid combination of clique and star expansions, although this involves the creation of additional nodes, and there may be scalability issues. In terms of predictive models, previous spectral-based hypergraph neural networks are analogous to applying GNNs on clique expansions, including HGNN (Feng et al., 2019), HCHA (Bai et al., 2021), H-GNNs (Zhang et al., 2022). Meanwhile, FastHyperGCN (Yadati et al., 2019) and HyperGCN (Yadati et al., 2019) reduce a hyperedge into a subgraph using Laplacian operators (Chan and Liang, 2020), which can be viewed as a modified form of clique expansion. HGAT (Ding et al., 2020), HNHN (Dong et al., 2020), HyperSAGE (Arya et al., 2020), UniGNN (Huang and Yang, 2021), (Srinivasan et al., 2021), Set-based models (Chien et al., 2022), (Heydari and Livi, 2022), Aponte et al. (2022), HEAT (Georgiev et al., 2022) take into account hyperedge features and use a message-passing framework, which can be interpreted as GNNs applied to the star expansion graph. And finally, (Wang et al., 2023) use gradient diffusion processes to motivate a broad class of hypergraph neural networks, although in the end there is not actually any specific energy function that is being minimized by the proposed model layers.

**Graph Neural Networks from Unfolded Optimization**. A variety of recent work has demonstrated that robust GNN architectures can be formed via graph propagation layers that mirror the unfolded descent iterations of a graph-regularized energy function (Chen and Eldar, 2021; Liu et al., 2021; Ma et al., 2020; Pan et al., 2021; Yang et al., 2021; Zhang et al., 2020; Zhu et al., 2021; Ahn et al., 2022). In doing so, the node embeddings at each layer can be viewed as increasingly refined approximations of an interpretable energy minimizer, that may be designed, for example, to mitigate GNN oversmoothing or perhaps inject robustness to spurious edges. Furthermore, these learnable embeddings can be integrated within a bilevel optimization framework (Wang et al., 2016) for supervised training. While at a high level we adopt a similar conceptual starting point, we nonetheless introduce non-trivial adaptations that are particular to the hypergraph domain, where this framework has not yet been extensively explored, and provide hypergraph-specific insights along the way.

## 3 Hypergraph Background and Notation

A hypergraph can be viewed as a higher-order form of graph whereby edges can encompass more than two nodes. Specifically, let \(\mathcal{G}(\mathcal{V},\mathcal{E})\) denote a hypergraph, where \(\mathcal{V}\) is a set of \(n=|\mathcal{V}|\) vertices and \(\mathcal{E}\) is a set of \(m=|\mathcal{E}|\) hyperedges. In contrast to a traditional graph, each hyperedge \(e_{k}\in\mathcal{E}\), can link an arbitrary number of nodes. The corresponding hypergraph connectivity structure is conveniently represented in a binary incidence matrix \(B\in\mathbb{R}^{n\times m}\), where \(B_{ik}=1\) if node \(v_{i}\in e_{k}\), otherwise \(B_{ik}=0\). We also use \(D_{H}\in\mathbb{R}^{m\times m}\) to denote the degree matrix of the hypergraph, where \(m_{e_{k}}\triangleq D_{H}[k,k]=\sum_{i}B_{ik}\).

And finally, we define input features and embeddings for both nodes and hyperedges. In this regard, \(X\in\mathbb{R}^{n\times d_{x}}\) represents a matrix of \(d_{x}\)-dimensional initial/given node features, while \(Y\in\mathbb{R}^{n\times d_{y}}\) refers to the corresponding node embeddings of size \(d_{y}\) we seek to learn. Analogously, \(U\in\mathbb{R}^{n\times d_{u}}\) and \(Z\in\mathbb{R}^{m\times d_{x}}\) are the initial edge features and learnable embeddings respectively. While here we have presented the most general form, we henceforth just assume \(d=d_{x}=d_{y}=d_{z}=d_{u}\) for simplicity.

## 4 A Family of Hypergraph Energy Functions

Our goal is to pursue hypergraph-based energy functions whose minima produce embeddings that will ultimately be useful for downstream predictive tasks. In this section, we first present an initial design of these functions followed by adaptations for handling the situation where no edge features \(U\) are available. We then show how in certain circumstances the proposed energy functions reduce to special cases that align with hypergraph star and clique expansions, before concluding with revised, simplified energy expressions informed by these considerations.

### Initial Energy Function Design and Motivation

We begin with the general form

\[\ell(Y,Z;\psi)=g_{1}(Y,X;\psi)+g_{2}(Z,U;\psi)+g_{3}(Y,Z,\mathcal{G};\psi) \tag{1}\]

where \(g_{1}(Y,X;\psi)\) and \(g_{2}(Z,U;\psi)\) are non-structural regularization factors over node and edge representations respectively, while \(g_{3}(Y,Z,\mathcal{G};\psi)\) explicitly incorporates hypergraph structure. In all cases \(\psi\) represents parameters that
control the shape of the energy, with particular choices that should be clear from the context (note that these parameters need not all be shared across terms; however, we nonetheless lump them together for notational convenience).

For the non-structural terms in (1), a plausible design criteria is to adopt functions that favor embeddings (either node or edge) that are similar to the corresponding input features or some transformation thereof. Hence we select

\[g_{1}(Y,X;\psi) = \sum_{i=1}^{n}\|y_{i}-f(x_{i};W_{x})\|_{2}^{2}\] \[g_{2}(Z,U;\psi) = \sum_{k=1}^{m}\|z_{k}-f(u_{k};W_{u})\|_{2}^{2}, \tag{2}\]

noting that both cases favor embeddings with minimal \(\ell_{2}\) distance from the trainable base predictor, and by extension, the initial features \(\{X,U\}\). In practice, the function \(f\) can be implemented as an MLP with node/edge weights \(W_{x}\) and \(W_{u}\) respectively.

Turning to \(g_{3}(Y,Z,\mathcal{G};\psi)\), our design is guided by the notion that:

1. Both node and edge embeddings should be individually constrained to a shared subset of \(\mathbb{R}^{d}\), e.g., consistent with most GNN architectures we may enforce non-negative embeddings;
2. Nodes sharing an edge should be similar when projected into an appropriate space, and;
3. Nodes within an edge set should have similar embeddings to the edge embedding, again, when suitably projected.

With these desiderata in mind, we adopt

\[g_{3}(Y,Z,\mathcal{G};\psi)=\sum_{i=1}^{n}\phi(y_{i})+\sum_{k=1} ^{m}\phi(z_{i})+\] \[\lambda_{0}\overbrace{\sum_{e_{k}\in\mathcal{E}}\sum_{i\in e_{k} }\sum_{j\in e_{k}}\sum_{j\in e_{k}}(|y_{i}H_{0}-y_{j}||_{2}^{2}}+\lambda_{1} \overbrace{\sum_{e_{k}\in\mathcal{E}}\sum_{i\in e_{k}}\sum_{i\in e_{k}}(|y_{i }H_{1}-z_{k}||_{2}^{2}}^{(b)} \tag{3}\]

For the first terms we choose \(\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}_{+}\) defined as \(\phi(p)\triangleq\sum_{j=1}^{d}\mathcal{I}_{\infty}[p_{i}<0]\), where \(\mathcal{I}_{\infty}\) is an indicator function that assigns an infinite penalty to any \(p_{i}<0\). This ensures that all node and edge embedding must be non-negative to achieve finite energy. Next, the term labeled (\(a\)) in (3) directly addresses criteria (ii). We note that the summation is over both indices \(i\) and \(j\) so that the symmetric counterpart, where the roles of nodes \(v_{i}\) and \(v_{j}\) are switched, is effectively included in the summation. And finally, criteria (iii) is handled by the last term, labeled (\(b\)). Here the node and edge embeddings play different roles and exhibit a natural asymmetry.1 Incidentally, the projections \(H_{0}\) and \(H_{1}\) can be viewed as compatibility matrices, initially introduced for label or belief propagation (Eswaran et al., 2017; Yamaguchi et al., 2016; Zhou et al., 2003) to provide additional flexibility to the metric in which entities are compared; for term (\(a\)) \(H_{0}\) facilitates the handling of nodes with potentially heterophily relationships, while for term (\(b\)) \(H_{1}\) accommodates the comparison of fundamentally different embedding types.

Footnote 1: While we could consider adding an additional factor \(||y_{i}-z_{k}H_{2}||_{2}^{2}\) to this term, we found that in practice it was not necessary.

### Handling a Lack of Edge Features

In some practical situations there may not be any initial hyperedge features \(U\). In such cases we could potentially modify \(\ell(Y,Z;\psi)\) accordingly in multiple different ways. First, and perhaps simplest, we can simply remove \(g_{2}(Z,U;\psi)\) from (1). We will explore the consequences of this option further in Section 4.3. But for tasks more related to hyperedge classification, it may be desirable to maintain this term for additional flexibility. Hence as a second option, we could instead create pseudo features \(\widetilde{U}\) with \(\widetilde{u}_{k}=\text{AGG}\left[\{x_{i}|i\in e_{k}\}\right]\) for all \(e_{k}\in\mathcal{E}\) for some aggregation function AGG. Or in a similar spirit, we could adopt \(f(u_{k};W_{u})\equiv\text{AGG}\left[\{f(x_{i};W_{x})|i\in e_{k}\}\right]\) such that aggregation now takes place after the initial feature transformations.

### Analysis of Simplified Special Cases

Because most hypergraph benchmarks for node classification, and many real-world use cases, involve data devoid of hyperedge features, in this section we more closely examine simplifications of (1) that arise when \(g_{2}(Z,U;\psi)\) is removed. For analysis purposes, it is useful to first introduce two representative hypergraph expansions, both of which can be viewed as converting the original hypergraph to a regular graph, which is tantamount to the assumption that edges in these expanded graphs involve only pairs of nodes.

**Clique Expansion.** For the _clique expansion_(Zien et al., 1999), we form the regular graph \(\mathcal{G}_{C}(\mathcal{V},\mathcal{E}_{C})\), where the node set \(\mathcal{V}\) remains unchanged while the edge set \(\mathcal{E}_{C}\) is such that, for all \(e_{k}\in\mathcal{E}\), we have that \(\{v_{i}|i\in e_{k}\}\) forms a complete subgraph of \(\mathcal{G}_{C}\). We define \(L_{C}\), \(A_{C}\), and \(D_{C}\) as the corresponding Laplacian, adjacency matrix, and degree matrix of \(\mathcal{G}_{C}\) respectively.

**Star Expansion.** In contrast, the _star expansion_(Zien et al., 1999) involves creating the bipartite graph \(\mathcal{G}_{S}(\mathcal{V}_{S},\mathcal{E}_{S})\), with revised node set \(\mathcal{V}_{S}=\{v_{1},\ldots,v_{n+m}\}\) and edge set \(\mathcal{E}_{S}\) defined such that \(\{v_{i},v_{n+k}\}\in\mathcal{E}_{S}\) iff
\(B_{ik}=1\). Conceptually, the resulting graph is formed with a new node associated with each hyperedge (from the original hypergraph), and an edge connecting every such new node to the original nodes within the corresponding hyperedges. Additionally, \(L_{S}=D_{S}-A_{S}\) is the revised Laplacian matrix, with \(D_{S}\) and \(A_{S}\) the degree and adjacent matrices of the star expansion graph.

**Unification.** We now introduce simplifying assumptions to link the proposed energy with the Laplacians of clique and star expansions as follows:

**Proposition 4.1**.: _Suppose \(g_{2}(Z,U;\psi)\) is removed from (1), \(H_{0}=H_{1}=I\), and define \(Z^{*}\triangleq D_{H}^{-T}B^{T}Y\). It then follows that_

\[\min_{Z}\ \ell(Y,Z;\psi) = g_{1}(Y,X;\psi)+\sum_{i=1}^{n}\phi(y_{i})\] \[+ 2\lambda_{0}\text{tr}[Y^{T}L_{C}Y]+\lambda_{1}\text{tr}\left( \left[\begin{array}{c}Y\\ Z^{*}\end{array}\right]^{T}L_{S}\left[\begin{array}{c}Y\\ Z^{*}\end{array}\right]\right)\] \[= g_{1}(Y,X;\psi)+\sum_{i=1}^{n}\phi(y_{i})\] \[+ 2\lambda_{0}\text{tr}[Y^{T}L_{C}Y]+\lambda_{1}\text{tr}[Y^{T} \bar{L}_{S}Y],\]

_where \(\bar{L}_{S}\triangleq\bar{D}_{S}-\bar{A}_{S}\), with \(\bar{A}_{S}\triangleq BD_{H}^{-1}B^{T}\) and \(\bar{D}_{S}\) a diagonal matrix with nonzero elements formed as the corresponding row-sums of \(\bar{A}_{S}\). Moreover, if \(\mathcal{G}\) is \(m_{e}\)-uniform,2 then under the same assumptions_

Footnote 2: An \(m_{e}\)-uniform hypergraph is such that every hyperedge joins exactly \(m_{e}\) nodes. Hence a regular graph is by default a 2-uniform hypergraph.

\[\min_{Z}\ \ell(Y,Z;\psi)=g_{1}(Y,X;\psi)+\sum_{i=1}^{n}\phi(y_{i})+\beta \text{tr}[Y^{T}L_{C}Y], \tag{5}\]

_where \(\beta\triangleq 2\lambda_{0}+\frac{\lambda_{1}}{m_{e}}\)._

All proofs are deferred to Appendix C. This last result demonstrates that, under the stated assumptions, the graph-dependent portion of the original hypergraph energy, after optimizing away the influence of \(Z\), can be reduced to a weighted quadratic penalty involving the graph Laplacian of the clique expansion. Moreover, this factor further resolves as

\[\text{tr}[Y^{T}L_{C}Y]=\frac{1}{2}\sum_{e_{k}\in\mathcal{E}}\sum_{i\in e_{k}} \sum_{j\in e_{k}}||y_{i}-y_{j}||_{2}^{2}. \tag{6}\]

Of course in more general settings, for example when \(H_{0}\neq H_{1}\neq I\), or when \(\phi(p)\neq\sum_{j=1}^{d}\mathcal{I}_{\infty}[p_{i}<0]\), this equivalence will _not_ generally hold.

### Revised Hypergraph Energy Functions

The analysis from the previous sections motivates two practical, revised forms of our original energy from (1), which we will later use for all of our empirical evaluations. For convenience, we define

\[\ell(Y;\psi)\triangleq\ell(Y,Z=Z^{*};\psi). \tag{7}\]

Then the first, more general variant, we adopt is

\[\ell(Y;\psi=\{W,H_{0},H_{1}\})\] \[= ||Y-f(X;W)||_{\mathcal{F}}^{2}+\sum_{i}\phi(y_{i})+\] \[\lambda_{0}\overbrace{\text{tr}\left[(YH_{0})^{T}D_{C}YH_{0}-2(YH _{0})^{T}A_{C}Y+Y^{T}D_{C}Y\right]}+\] \[\lambda_{1}\overbrace{\text{tr}\left[(YH_{1})^{T}\bar{D}_{S}YH_{ 1}-2(YH_{1})^{T}BZ^{*}+Z^{*T}D_{H}Z^{*}\right]},\]

where \(\bar{D}_{S}\) is defined as in Proposition 4.1. Moreover, to ease later exposition, we have overloaded the definition of \(f\) such that \(\|Y-f(X;W)\|_{\mathcal{F}}^{2}\equiv\sum_{i=1}^{n}\|y_{i}-f(x_{i};W)\|_{2}^{2}\). And secondly, as a less complex alternative we have

\[\ell(Y;\psi=\{W,I,I\})= \tag{9}\] \[||Y-f(X;W)||_{\mathcal{F}}^{2}+\sum_{i}\phi(y_{i})+\text{tr}[Y^{T }(\lambda_{0}L_{C}+\lambda_{1}\bar{L}_{S})Y].\]

## 5 Hypergraph Node Classification via Bilevel Optimization

We now demonstrate how the optimal embeddings obtained by minimizing the energy functions from the previous section can be applied to our ultimate goal of hypergraph node classification. For this purpose, define

\[Y^{*}(\psi)=\arg\min_{Y}\ell(Y;\psi), \tag{10}\]

noting that the solution depends explicitly on the parameters \(\psi\) governing the shape of the energy. We may then consider treating \(Y^{*}(\psi)\), which is obtainable from the above optimization process, as features to be applied to a discriminative node classification loss \(\mathcal{D}\) that can be subsequently minimized via a second, meta-level optimization step.3 In aggregate we arrive at the _bilevel_ optimization problem

Footnote 3: Because our emphasis is hypergraph node classification, we will not explicitly use any analogous hyperedge embeddings for the meta-level optimization; however, they nonetheless still play a vital role given that they are co-adapted with the node embeddings during the lower-level optimization per the discussion from the previous section.

\[\ell(\theta,\psi)\triangleq\sum_{i=1}^{n^{\prime}}\mathcal{D}(h[y_{i}^{*}( \psi);\theta],\tau_{i}), \tag{11}\]
where \(\mathcal{D}\) is chosen as an classification-friendly cross-entropy function, \(y_{i}^{*}(\psi)\) is the \(i\)-th row of \(Y^{*}(\psi)\), and \(\tau_{i}\in\mathbb{R}^{c}\) is the ground-truth label of node \(i\) to be approximated by some differentiable node-wise function \(h:\mathbb{R}^{d}\rightarrow\mathbb{R}^{c}\) with trainable parameters \(\theta\). We have also implicitly assumed that the first \(n^{\prime}\) nodes of \(\mathcal{G}\) are labeled. Intuitively, (11) involves training a classifier \(h\), with input features \(y_{i}^{*}(\psi)\), to predict labels \(\tau_{i}\).

At this point, assuming \(\partial Y^{*}(\psi)/\partial\psi\) is somehow computable, then \(\ell(\psi,\theta)\) can be efficiently trained over _all_ parameters, including \(\psi\) from the lower level optimization. However, directly computing \(\partial Y^{*}(\psi)/\partial\psi\) is not generally feasible. Instead, in the remainder of this section we will derive approximate embeddings \(\hat{Y}(\psi)\approx Y^{*}(\psi)\) whereby \(\partial\hat{Y}(\psi)/\partial\psi\) can be computed efficiently. And as will be assessed in greater detail later, the computational steps we derive to produce \(\hat{Y}(\psi)\) will mirror the layers of canonical graph neural network architectures. It is because of this association that we refer to our overall model as **PhenomNN**, for _Purposeful Hyper-Edges iN Optimization Motivated Neural Networks_ as mentioned in the introduction.

### Deriving Proximal Gradient Descent Steps

To efficiently deploy proximal gradient descent (PGD) (Parikh et al., 2014), we first must split our loss into a smooth, differentiable part, and a non-smooth but separable part. Hence we adopt the decomposition

\[\ell(Y;\psi)=\bar{\ell}(Y;\psi)+\sum_{i}\phi(y_{i}), \tag{12}\]

where \(\bar{\ell}(Y;\psi)\) is defined by exclusion upon examining the original form of \(\ell(Y;\psi)\). The relevant proximal operator is

\[\textbf{prox}_{\phi}(V) \triangleq \arg\min_{Y}\frac{1}{2}||V-Y||_{\mathcal{F}}^{2}+\sum_{i}\phi(y_ {i}) \tag{13}\] \[= \text{max}(0,V),\]

where the max operator is assumed to apply elementwise. Subsequent PGD iterations for minimizing (12) are then computed as

\[\tilde{Y}^{(t+1)} =Y^{(t)}-\alpha\Omega\nabla_{Y^{(t)}}\bar{\ell}(Y^{(t)};\psi) \tag{14}\] \[Y^{(t+1)} =\text{max}(0,\tilde{Y}^{(t+1)}), \tag{15}\]

where \(\alpha\) is a step-size parameter and \(\Omega\) is a positive-definite pre-conditioner to be defined later. Incidentally, as will become apparent shortly, (14) will occupy the role of a pre-activation hypergraph neural network layer, while (15) provides a ReLU nonlinearity. A related association was previously noted within the context of traditional GNNs (Yang et al., 2021). We now examine two different choices for \(\Omega\) and \(\psi\) that correspond with the general form from (8) and the simplified alternative from (9).

General FormTo compute (14), we consider term (\(a\)) and (\(b\)) from (8) separately. Beginning with (\(a\)), the corresponding gradient is

\[2D_{C}Y-2\tilde{Y}_{C}, \tag{16}\]

where \(\tilde{Y}_{C}\triangleq A_{C}Y(H_{0}+H_{0}^{T})-D_{C}YH_{0}H_{0}^{T}\). Similarly, for (\(b\)) the gradient is given by

\[2BD_{D}^{-1}D_{H}(BD_{H}^{-1})^{T}Y-2\tilde{Y}_{S}, \tag{17}\]

where \(\tilde{Y}_{S}\triangleq(B(BD_{H}^{-1})^{T}YH_{1}^{T}+BD_{H}^{-1}B^{T}YH_{1})- \bar{D}_{S}YH_{1}H_{1}^{T}\). Additionally, given that \(BD_{H}^{-1}D_{H}(BD_{H}^{-1})^{T}=B(BD_{H}^{-1})^{T}=BD_{H}^{-1}B^{T}=\bar{A} _{S}\), we can reduce (17) to

\[2\bar{A}_{S}Y-2\tilde{Y}_{S}, \tag{18}\]

since now \(\tilde{Y}_{S}=\bar{A}_{S}Y(H_{1}+H_{1}^{T})-\bar{D}_{S}YH_{1}H_{1}^{T}\). Combining terms, the gradient for \(\bar{\ell}(Y;\psi)\) is

\[\frac{\bar{\ell}(Y;\psi)}{\partial Y}=2\lambda_{0}(D_{C}Y-\tilde{ Y}_{C})+2\lambda_{1}(\bar{A}_{S}Y-\tilde{Y}_{S})\] \[+2Y-2f\left(X;W\right), \tag{19}\]

and (14) becomes

\[\bar{Y}^{(t+1)}=Y^{(t)}-\alpha\Bigg{[}\lambda_{0}(D_{C}Y^{(t)}- \tilde{Y}_{C}^{(t)})\] \[+\lambda_{1}(\bar{A}_{S}Y^{(t)}-\tilde{Y}_{S}^{(t)})+Y^{(t)}-f \left(X;W\right)\Bigg{]}, \tag{20}\]

where \(\alpha/2\) is the step size. The coefficient \(\bar{\Omega}\) before \(Y^{(t)}\) is

\[\bar{\Omega}\triangleq\lambda_{0}D_{C}+\lambda_{1}\bar{A}_{S}+I. \tag{21}\]

Applying Jacobi preconditioning (Axelsson, 1996) often aids convergence by helping to normalize the scales across different dimensions. One natural candidate for the preconditioner is \(\left(\text{diag}[\bar{\Omega}]\right)^{-1}\); however, we use the more spartan \(\Omega=\tilde{D}^{-1}\) where \(\tilde{D}\triangleq\lambda_{0}D_{C}+\lambda_{1}\bar{D}_{S}+I\). After rescaling and applying (15), the composite PhenomNN update is given by

\[Y^{(t+1)} =\text{ReLU}\Bigg{(}(1-\alpha)Y^{(t)}+\alpha\tilde{D}^{-1}\Big{[} f\left(X;W\right)\] \[+\lambda_{0}\tilde{Y}_{C}^{(t)}+\lambda_{1}(\bar{L}_{S}Y^{(t)}+ \tilde{Y}_{S}^{(t)})\Big{]}\Bigg{)}, \tag{22}\]

where \(\bar{L}_{S}=\bar{D}_{S}-\bar{A}_{S}\) as in Proposition 4.1. This respresents the general form of PhenomNN.

Simplified AlternativeRegarding the simplified energy from (9), the relevant gradient is

\[\frac{\partial\bar{\ell}(Y;\psi=W,I,I)}{\partial Y}=2(\lambda_{0} L_{C}+\lambda_{1}\bar{L}_{S})Y+\] \[2Y-2f\left(X;W\right), \tag{23}\]
leading to the revised update

\[\bar{Y}^{(t+1)}=Y^{(t)}-\alpha\left[\tilde{\Omega}Y^{(t)}-f\left(X;W \right)\right],\] \[\text{with}\;\;\tilde{\Omega}\triangleq\lambda_{0}L_{C}+\lambda_{ 1}\bar{L}_{S}+I \tag{24}\]

and step size \(\alpha/2\) as before. And again, we can apply preconditioning, in this case rescaling each gradient step by \(\Omega=\left(\text{diag}[\tilde{\Omega}]\right)^{-1}=\left(\lambda_{0}D_{C}+ \lambda_{1}\bar{D}_{S}+I\right)^{-1}=\tilde{D}^{-1}\). So the final/composite update formula, including (15), becomes

\[Y^{(t+1)}=\text{ReLU}\Big{(}(1-\alpha)Y^{(t)}+\] \[\alpha\tilde{D}^{-1}\left[(\lambda_{0}A_{C}+\lambda_{1}\bar{A}_{S })Y^{(t)}+f\left(X;W\right)\right]\Big{)}. \tag{25}\]

We henceforth refer to this variant as **PhenomNN\({}_{\text{simple}}\)**.

### Overall Algorithm

The overall algorithm for PhenomNN is demonstrated in Algorithm 1.

```
Input: Hypergraph incidence matrix \(B\), node features \(X\), number of layers \(T\), training epochs \(E\), and node labels \(\tau=\{\tau_{i}\}\). for\(e=0\) to \(E-1\)do  Set initial projection \(Y^{(0)}=f(X;W)\), where \(f\) is the trainable base model. for\(t=0\) to \(T-1\)do \(Y^{(t+1)}=Update(Y^{(t)})\), where \(Update\) is computed via (22) for PhenomNN or (25) for PhenomNN\({}_{\text{simple}}\). endfor  Compute loss \(\ell(\theta,\psi)=\sum_{i}\mathcal{D}(h[y_{i}^{(T)};\theta],\tau_{i})\) from (11), where \(\psi=\{W,H_{0},H_{1}\}\) for PhenomNN and \(\psi=\{W,I,I\}\) for PhenomNN\({}_{\text{simple}}\), noting that each \(y_{i}^{(T)}\) is a trainable function of \(\psi\) by design.  Backpropagate over all parameters \(\psi,\theta\) using optimizer (Adam, SGD, etc.) endfor
```

**Algorithm 1** PhenomNN Algorithm for Hypergraph Node Classification.

### Convergence Analysis

We now consider the convergence of the iterations (22) and (25) introduced in the previous section. First, for the more general form we have the following:

**Proposition 5.1**.: _The PhenomNN updates from (22) are guaranteed to monotonically converge to the unique global minimum of \(\ell(Y;\psi)\) on the condition that_

\[\alpha<\frac{1+\lambda_{0}d_{\text{Cmin}}+\lambda_{1}d_{\text{Smin}}}{1+ \lambda_{0}d_{\text{Cmin}}+\sigma_{\text{max}}}, \tag{26}\]

_where \(d_{\text{Cmin}}\) is the minimum diagonal element of \(I\otimes D_{C}\), \(d_{\text{Smin}}\) is the minimum diagonal element of \(I\otimes\bar{D}_{S}\) and \(\sigma_{\text{max}}\) is the max eigenvalue of \((Q-P+\lambda_{1}I\otimes\bar{A}_{S})\) with_

\[Q \triangleq\lambda_{0}H_{0}^{T}H_{0}\otimes D_{C}+\lambda_{1}H_{1}^ {T}H_{1}\otimes\bar{D}_{S}, \tag{27}\] \[P \triangleq\lambda_{0}(H_{0}+H_{0}^{T})\otimes A_{C}+\lambda_{1}(H _{1}+H_{1}^{T})\otimes\bar{A}_{S}. \tag{28}\]

And for the restricted case where \(\psi=\{W,I,I\}\), the convergence conditions simplify as follows:

**Corollary 5.2**.: _The PhenomNN\({}_{\text{simple}}\) updates from (25) are guaranteed to monotonically converge to the unique global minimum of \(\ell(Y;\psi=\{W,I,I\})\) on the condition that_

\[\alpha<\frac{1+\lambda_{0}d_{\text{Cmin}}+\lambda_{1}d_{\text{Smin}}}{1+ \lambda_{0}d_{\text{Cmin}}+\lambda_{1}d_{\text{Smin}}-\sigma_{\text{min}}}, \tag{29}\]

_where \(\sigma_{min}\) is the min eigenvalue of \((\lambda_{0}A_{C}+\lambda_{1}\bar{A}_{S})\)._

### Complexity Analysis

Analytically, PhenomNN\({}_{\text{simple}}\) has a time complexity given by \(O(|\mathcal{E}|Td+|\mathcal{V}|Pd^{2})\), where \(|\mathcal{E}|\) is edge number, \(|\mathcal{V}|\) is the node number, \(T\) is the number of layers/iterations, \(d\) is the hidden size, and \(P\) is the number of MLP layers in \(f(\cdot;W)\). In contrast, for PhenomNN this complexity increases to \(O(|\mathcal{E}|Td+|\mathcal{V}|(T+P)d^{2})\), which is roughly the same as a standard GCN model. In fact, the widely-used graph convolution networks (GCN) (Kipf and Welling, 2016) have equivalent complexity to PhenomNN up to the factor of \(P\) which is generally small (e.g., \(P=1\) for PhenomNN in our experiments, while for a GCN \(P=0\)). In this way then, PhenomNN\({}_{\text{simple}}\) is actually somewhat cheaper than a GCN when \(T>P\). Additionally, we include complementary empirical results related to time and space complexity in Section 7.

## 6 Connections with Existing GNN Layers

As mentioned in Section 4.3, the clique and star expansions can be invoked to transform hypergraphs into homogeneous and bipartite graphs respectively (where the latter is a special case of a heterogeneous graph). In this section we examine how the layer-wise structure of two of the most popular GNN models, namely GCN (Kipf and Welling, 2016) mentioned previously, and relational graph convolution networks (RGCN) (Schlichtkrull et al., 2018), relate to PhenomNN and simplifications thereof.

### Homogeneous Graphs and GCN

Using the so-called message-passing form of expression, the embedding update for the \(i\)-th node of the \(t\)-th GCN
layer can be written as

\[y_{i}^{(t+1)}=\sigma\left(\sum_{j\in\mathcal{N}_{i}}\frac{1}{c_{ij}}W^{(t)}y_{j}^{ (t)}\right) \tag{30}\]

where \(\sigma\) is an activation function like ReLU, \(W^{(t)}\) are weights, \(c_{ij}\triangleq\sqrt{|\mathcal{N}_{i}||\mathcal{N}_{j}|}\) and \(\mathcal{N}_{i}\) refers to the set of neighboring nodes in some input graph (note also that the graph could have self-loops in which case \(i\in\mathcal{N}_{i}\)). Interestingly, follow-up work (Ma et al., 2020; Pan et al., 2021; Yang et al., 2021; Zhang et al., 2020; Zhu et al., 2021) has demonstrated that this same basic layer-wise structure can be closely linked to iterative steps designed to minimize the energy

\[\ell(Y)=||Y-f(X;W)||_{\mathcal{F}}^{2}+\lambda\text{tr}[Y^{T}LY], \tag{31}\]

where \(f\) is defined as before and \(L\) is the assumed graph Laplacian matrix. One way to see this is to examine a preconditioned gradient step along (31), which can be expressed as

\[Y^{(t+1)}=(1-\alpha)Y^{(t)}+\alpha\widetilde{D}_{0}^{-1}\left[ \lambda AY^{(t)}+f\left(X;W\right)\right], \tag{32}\]

with preconditioner \(\widetilde{D}_{0}^{-1}=(\lambda D+I)^{-1}\), step-size parameter \(\alpha\), graph adjacency matrix \(A\), and corresponding degree matrix \(D\). Moreover, for a single node \(i\), (32) can be reduced to

\[y_{i}^{(t+1)}=\left(\sum_{j\in\mathcal{N}_{i}}\frac{1}{\tilde{c}_{i}}y_{j}^{( t)}\right)+\tilde{f}_{i}(x_{i};W), \tag{33}\]

where \(\tilde{c}_{i}\) is a scaling constant dependent on \(\lambda\), the gradient step-size, and the preconditioner, while \(\tilde{f}_{i}\) is merely \(f\) similarly rescaled. If we add an additional penalty \(\phi\) and subsequent proximal operator step to introduce a non-linearity, then this result is very similar to (30), although without the weight matrix directly on each \(y_{j}^{(t)}\) but with an added skip connection to the input layer.

Importantly for our purposes though, if the input graph is chosen to be a hypergraph clique expansion, and we set \(D=D_{C}\), \(A=A_{C}\), \(\lambda=\lambda_{0}\), and \(\lambda_{1}=0\), then we arrive at a special case of PhenomNNsimple from (25). Of course one might not naturally conceive of the more generalized form that leads to PhenomNNsimple, and by extension PhenomNN, without the interpretable grounding of the underlying hypergraph energy functions involved.

### Heterogeneous Graphs and RGCN

For heterogeneous graphs applied to RGCN, the analogous message-passing update for the \(i\)-th node in the \(t\)-th layer is given by

\[y_{i}^{(t+1)}=\sigma\left(\sum_{r\in\mathcal{R}}\sum_{j\in\mathcal{N}_{i}^{r} }\frac{1}{c_{i,r}}y_{j}^{(t)}W_{r}^{(t)}+y_{i}^{(t)}W_{0}^{(t)}\right), \tag{34}\]

where \(\mathcal{R}\) is the set of edge types in a heterogeneous input graph, \(\mathcal{N}_{i}^{r}\) is the set of neighbors with edge type \(r\), \(c_{i,r}\triangleq|\mathcal{N}_{i}^{r}|\), and \(W_{r}^{(t)}\) and \(W_{0}^{(t)}\) are weight/projection matrices. In this context, the RGCN input could conceivably be chosen as the bipartite graph produced by a given star expansion (e.g., such a graph could be assigned the edge types "hypergraph node belongs to hyperedge" and "hyper-edge belongs to hypergraph node").

For comparison purposes, we can also re-express our general PhenomNN model from (22), in the node-wise message-passing form

\[y_{i}^{(t+1)}=\sigma\left(\sum_{j\in\mathcal{N}_{i}^{C}}y_{j}^{( t)}W_{ij}^{(t)}+y_{i}^{(t)}W_{i}^{(t)}+\alpha\tilde{D}_{ii}^{-1}f\left(x_{i};W \right)\right), \tag{35}\]

where \(\mathcal{N}_{i}^{C}\) are neighbors in the clique (not star) expansion graph (more on this below) and the weight matrices are characterized by the special energy-function-dependent forms

\[W_{ij}^{(t)} \triangleq\alpha\tilde{D}_{ii}^{-1}\Big{[}\lambda_{0}A_{C}[i,j](H_ {0}+H_{0}^{T})+\] \[\lambda_{1}\bar{A}_{S}[i,j](H_{1}+H_{1}^{T}-I)\Big{]}, \tag{36}\] \[W_{i}^{(t)} \triangleq(1-\alpha)I-\alpha\tilde{D}_{ii}^{-1}\Big{[}\lambda_{0} D_{C}[i,i]H_{0}H_{0}^{T}+\] \[\lambda_{1}\bar{D}_{S}[i,i](H_{1}H_{1}^{T}-I)\Big{]}. \tag{37}\]

While the basic structures of (34) and (35) are similar, there are several key differences:

* When RGCN is applied to the star expansion, neighbors are defined by the resulting bipartite graph, and nodes in the original hypergraph do not directly pass messages to each other. In contrast, because within PhenomNN we have optimized away the hyperedge embeddings, the implicit graph that dictates neighborhood structure is actually the _clique expansion graph_ as reflected in (35).
* The PhenomNN projection matrices have special structure infused from the energy function and optimization over the edge embeddings. As such, unlike RGCN node \(i\) receives messages from its connected neighbors and itself, with projection matrices \(W_{ij}^{(t)}\) and \(W_{i}^{(t)}\) that can vary from node to node and edge to edge. In contrast, RGCN has layer-wise (or analogously iteration-wise) dependent weights.
* PhenomNN has an additional weighted skip connection from the input base model \(f\left(x_{i};W\right)\). While of course RGCN could also be equipped with a similar term, this would be accomplished in a post-hoc fashion, and not tethered to an underlying energy function.

## 7 Hypernode Classification Experiments

In this section we evaluate PhenomNNsimple and PhenomNN on various hypergraph benchmarks focusing on hypernode classification and compare against previous SOTA approaches.

Datasets.Existing hypergraph benchmarks mainly focus on hypernode classification. We adopt five public citation network datasets from (Zhang et al., 2022): Co-authorship/Cora,Co-authorship/DBLP, Co-citaion/Cora, Co-citaion/Pubmed, Co-citaion/Citeseer. These datasets and splits are constructed by (Yadati et al., 2019) ([https://github.com/malllabiisc/HyperGCN](https://github.com/malllabiisc/HyperGCN)). We also adopt two other public visual object classification datasets: Princeton ModelNet40 (Wu et al., 2015) and the National Taiwan University (NTU) 3D model dataset (Chen et al., 2003). We follow HGNN (Feng et al., 2019) to preprocess the data by MVCNN (Su et al., 2015) and GVCNN (Feng et al., 2018) and obtain the hypergraphs. Additionally, we use the datasets provided by the public code ([https://github.com/iMoonLab/HGNN](https://github.com/iMoonLab/HGNN)) associated with (Feng et al., 2019). Finally, (Chien et al., 2022) construct a public hypergraph benchmark for hypernode classification which includes ModelNet40\({}^{*}\), NTU2012\({}^{*}\), Yelp (Yelp), House (Chodrow et al., 2021), Walmart (Amburg et al., 2020), and 20News (Dua and Graff, 2017). ModelNet40\({}^{*}\) and NTU2012\({}^{*}\) have the same raw data as ModelNet40 and NTU2012 mentioned before in (Zhang et al., 2022) but different splits. All datasets from (Chien et al., 2022) are downloaded from their code site ([https://github.com/jianhao2016/AllSet](https://github.com/jianhao2016/AllSet)).4

Footnote 4: Note that we excluded a few datasets for the following reasons: The Zoo dataset is very small; the Mushroom dataset is too easy; the Citation datasets are similar to (Zhang et al., 2022), and since we have ModelNet40* and NTU2012\({}^{*}\) for comparison of different baselines from both papers, we did not select them.

Baselines.For datasets from (Zhang et al., 2022), we adopt the baselines from their paper which includes a multi-layer perceptron with explicit hypergraph Laplacian regularization (MLP+HLR), FastHyperGCN (Yadati et al., 2019), HyperGCN (Yadati et al., 2019), HGNN (Feng et al., 2019), HNHN (Dong et al., 2020), HGAT (Ding et al., 2020), HyperSAGE (Arya et al., 2020), UniGNN (Huang and Yang, 2021), and various hypergraph GNNs (H-GNNs) (Zhang et al., 2022) proposed by them. For datasets from (Chien et al., 2022), we also select baselines from their paper including an MLP, CE (Clique Expansion\(+\)GCN, CE\(+\)GAT, HNHN, HGNN, HCHA (Bai et al., 2021), HyperGCN, UniGCNII (Huang and Yang, 2021), HAN (Wang et al., 2019) with full batch and mini-batch settings, and AllsetTransformer and AllDeepSets (Chien et al., 2022).

Implementations.We use a one-layer MLP for \(f(X;W)\). Also, in practice we found that only using ReLU at the end of propagation steps works well. Detailed hyperparameter settings are deferred to Appendix D. We choose the hidden dimension of our models to be the same or less than the baselines in previous work. For results in Table 1, we conduct experiments on 10 different train-test splits and report average accuracy of test samples following (Zhang et al., 2022). For results in Table 2, we randomly split the data into training/validation/test samples using (50%/25%/25%) splitting percentages as in (Chien et al., 2022) and report

\begin{table}
\begin{tabular}{c|c c c c c c c c} \hline \hline  & Cora & DBLP & Cora & Pubmed & Citeseer & NTU2012 & ModelNet40 & Avg Ranking \\  & (co-authorship) & (co-authorship) & (co-citation) & (co-citation) & (co-citation) & (both features) & (both features) & Avg Ranking \\ \hline MLP+HLR & 59.8 \(\pm\) 4.7 & 63.6 \(\pm\) 7.4 & 61.0 \(\pm\) 4.1 & 64.7 \(\pm\) 3.1 & 56.4 \(\pm\) 2.6 & - & - & 13.6 \\ FastHyperGCN & 61.1 \(\pm\) 8.2 & 68.1 \(\pm\) 9.6 & 61.3 \(\pm\) 10.3 & 65.7 \(\pm\) 11.1 & 56.2 \(\pm\) 8.1 & - & - & 12.4 \\ HyperGCN & 63.9 \(\pm\) 7.3 & 70.9 \(\pm\) 8.3 & 62.5 \(\pm\) 9.7 & 67.8 \(\pm\) 9.5 & 57.3 \(\pm\) 7.3 & - & - & 10.8 \\ HGNN & 63.2 \(\pm\) 3.1 & 68.1 \(\pm\) 9.6 & 70.9 \(\pm\) 2.9 & 68.7 \(\pm\) 8.7 & 56.7 \(\pm\) 3.8 & 83.54 \(\pm\) 0.50 & 97.15 \(\pm\) 0.14 & 9.4 \\ HBNN & 64.0 \(\pm\) 2.4 & 84.4 \(\pm\) 3.0 & 41.6 \(\pm\) 3.1 & 41.9 \(\pm\) 4.7 & 33.6 \(\pm\) 2.1 & - & - & 13.0 \\ HGAT & 65.4 \(\pm\) 1.5 & OOM & 52.2 \(\pm\) 3.5 & 46.3 \(\pm\) 0.5 & 38.3 \(\pm\) 1.5 & 84.05 \(\pm\) 0.36 & 96.44 \(\pm\) 0.15 & 12.0 \\ HyperSAGE & 72.4 \(\pm\) 1.6 & 77.4 \(\pm\) 3.8 & 69.3 \(\pm\) 2.7 & 72.9 \(\pm\) 1.3 & 61.8 \(\pm\) 2.3 & - & - & 8.6 \\ UniGNN & 75.3 \(\pm\) 1.2 & 88.8 \(\pm\) 0.2 & 70.1 \(\pm\) 1.4 & 74.4 \(\pm\) 1.0 & 63.6 \(\pm\) 1.3 & 84.55 \(\pm\) 0.40 & 96.69 \(\pm\) 0.07 & 6.0 \\ \hline H-ChebNet & 70.6 \(\pm\) 2.1 & 87.9 \(\pm\) 0.24 & 69.7 \(\pm\) 2.0 & 74.3 \(\pm\) 1.5 & 63.5 \(\pm\) 1.3 & 83.16 \(\pm\) 0.46 & 96.95 \(\pm\) 0.09 & 8.0 \\ H-APPNP & 76.4 \(\pm\) 0.8 & 89.4 \(\pm\) 0.18 & 70.9 \(\pm\) 0.7 & 75.3 \(\pm\) 1.1 & 64.5 \(\pm\) 1.4 & 83.57 \(\pm\) 0.42 & 97.20 \(\pm\) 0.14 & 4.6 \\ H-SSFC & 72.0 \(\pm\) 1.2 & 88.6 \(\pm\) 0.16 & 68.8 \(\pm\) 2.1 & 74.5 \(\pm\) 1.3 & 60.5 \(\pm\) 1.7 & 84.13 \(\pm\) 0.34 & 97.07 \(\pm\) 0.07 & 7.6 \\ H-GCN & 74.8 \(\pm\) 0.9 & 89.0 \(\pm\) 0.19 & 69.5 \(\pm\) 2.0 & 75.4 \(\pm\) 1.2 & 62.7 \(\pm\) 1.2 & 84.5 \(\pm\) 0.40 & 97.28 \(\pm\) 0.15 & 5.4 \\ H-GCNII & 76.2 \(\pm\) 1.0 & 89.8 \(\pm\) 0.20 & 72.5 \(\pm\) 1.2 & 75.8 \(\pm\) 1.1 & 64.5 \(\pm\) 1.0 & 85.17 \(\pm\) 0.36 & 97.5 \(\pm\) 0.07 & 3.0 \\ \hline PhenomNNsimple & **77.62 \(\pm\) 1.30** & 89.74 \(\pm\) 0.16 & 72.81 \(\pm\) 1.67 & 76.20 \(\pm\) 1.41 & 65.07 \(\pm\) 1.08 & 85.39 \(\pm\) 0.40 & **97.83 \(\pm\) 0.09** & 1.9 \\ PhenomNN & **77.11 \(\pm\) 0.45** & **98.31 \(\pm\) 0.05** & **73.09 \(\pm\) 0.05** & **78.12 \(\pm\) 0.24** & **65.77 \(\pm\) 0.45** & **85.40 \(\pm\) 0.42** & 97.77 \(\pm\) 0.11 & **13** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Results on datasets from (Zhang et al., 2022): Mean accuracy (%) ± standard deviation results over 10 train-test splits. Boldfacters are used to indicate the best mean accuracy and underline for the second. “-” means not reported in their paper so in average ranking we just average over the ones that are available. OOM indicates out-of-memory.
the average accuracy over ten random splits. All experiments are implemented on RTX 3090 with Pytorch and DGL (Wang et al., 2019).

**Results.** As shown in Table 1, our models achieve the best performance and top ranking on all datasets from (Zhang et al., 2022) compared to previous baselines. And in Table 2, our models achieve the first (PhenomNN\({}_{\text{simple}}\)) and tied-for-second (PhenomNN) overall performance ranking on the benchmarks from (Chien et al., 2022).

**Empirical evaluation of time and space complexity.** In practice, we find that PhenomNN is roughly \(2\times\) to \(3\times\) slower than a GCN given the integration of two expansions based on \(H_{0}\) and \(H_{1}\), which implies that the constant multiplying the theoretical complexity from above is at least doubled as expected. Of course timing results will still vary based on hardware and implementation details. As an example, we measure the training time of GCN and our models on the same hardware on Coauthorship-DBLP data with hidden size 64 and 8 layers. We observe 0.047s/epoch for GCN and 0.045s/epoch for PhenomNN\({}_{\text{simple}}\) and 0.143s/epoch for PhenomNN under these conditions. In terms of the space efficiency, our models are also analytically similar to common GNNs. And under the same settings as above, the memory consumption is 1665MB for GCN, 1895MB for PhenomNN\({}_{\text{simple}}\), and 2424MB for PhenomNN.

**Ablations.** For space considerations, we defer ablations to Appendix B; however, we nonetheless highlight some of our findings here. For example, in Table 5 (Appendix B) we demonstrate the effect of different hypergraph energy function terms, which are associated with different hypergraph expansions per Proposition 4.1. In brief here, we explore different selections of \(\{\lambda_{0},\lambda_{1}\}\in\{\{0,1\},\{1,0\},\{1,1\},\}\) which in effect modulate the inclusion of clique- and star-like expansion factors. Results demonstrate that on most datasets, the combination of both expansions, with their complementary roles, is beneficial.

We also explore the tolerance of our model to different hidden dimensions in Table 6. In brief, we fix other hyperparameters and obtain results across different hidden dimensions with PhenomNN\({}_{\text{simple}}\) for simplicity; results for PhenomNN are similar. Overall, this ablation demonstrates the stability of our approach across hidden dimension.

**Additional comparisons and discussion.** As suggested by reviewers, we include additional discussion and comparison with existing work in Appendix A due to the page limit. This includes side-by-side evaluations with RGCN and the model from (Wang et al., 2023) which was not yet published at the time of our original submission.

## 8 Conclusion

While hypergraphs introduce compelling modeling flexibility, they still remain relatively under-explored in the GNN literature. With the potential to better understand hypergraph properties and expand their utility, we have introduced an expressive family of hypergraph energy functions and fleshed out their connection with previous hypergraph expansions. We then leverage this perspective to design what can be interpreted as hypergraph neural network layers that are in one-to-one correspondence with proximal gradient steps descending these energies. We also characterize the similarities and differences of these layers w.r.t. popular existing GNN architectures. In the end, the proposed framework achieves competitive or SOTA performance on key hypergraph node classification benchmarks.

## 9 Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 62236004 and No. 62022027)

\begin{table}
\begin{tabular}{c|c c c c c c c c|c} \hline  & NTU2012\({}^{*}\) & ModelNet40\({}^{*}\) & Yelp & House(1) & Walmart(1) & House(0.6) & Walmart(0.6) & 20Newsgroups & Avg Ranking \\ \hline MLP & \(85.52\pm 1.49\) & \(96.14\pm 0.36\) & \(31.96\pm 0.44\) & \(67.93\pm 2.33\) & \(45.51\pm 0.24\) & \(81.53\pm 2.26\) & \(63.28\pm 0.37\) & \(81.42\pm 0.49\) & 6.9 \\ CEGCN & \(81.52\pm 1.43\) & \(89.92\pm 0.46\) & OOM & \(62.80\pm 2.61\) & \(54.44\pm 0.24\) & \(64.36\pm 2.41\) & \(59.78\pm 0.32\) & OOM & 11.5 \\ CEGAT & \(82.21\pm 1.23\) & \(92.52\pm 0.39\) & OOM & \(69.09\pm 0.39\) & \(51.14\pm 0.56\) & \(77.18\pm 2.53\) & \(59.47\pm 1.05\) & OOM & 10.4 \\ HNHN & \(89.11\pm 1.44\) & \(97.84\pm 0.25\) & \(31.65\pm 0.44\) & \(67.80\pm 2.59\) & \(47.18\pm 0.35\) & \(78.78\pm 1.88\) & \(65.80\pm 0.39\) & \(31.35\pm 0.61\) & 7.1 \\ HGNN & \(87.72\pm 1.35\) & \(95.44\pm 0.33\) & \(33.04\pm 0.62\) & \(61.39\pm 0.29\) & \(62.00\pm 0.24\) & \(66.16\pm 1.80\) & \(77.72\pm 0.21\) & \(80.33\pm 0.42\) & 7.8 \\ HCHA & \(87.48\pm 1.87\) & \(94.48\pm 0.28\) & \(30.99\pm 0.72\) & \(61.36\pm 2.53\) & \(62.45\pm 0.26\) & \(67.91\pm 2.26\) & \(77.12\pm 0.26\) & \(80.33\pm 0.80\) & 8.8 \\ HyperCCN & \(56.36\pm 4.86\) & \(75.58\pm 2.56\) & \(29.42\pm 1.54\) & \(48.31\pm 2.39\) & \(44.74\pm 2.81\) & \(78.22\pm 2.46\) & \(55.31\pm 0.30\) & \(81.05\pm 0.59\) & 12 \\ UniGCN & \(89.30\pm 1.33\) & \(89.07\pm 0.23\) & \(31.70\pm 0.52\) & \(67.25\pm 2.57\) & \(54.45\pm 0.37\) & \(68.05\pm 1.96\) & \(72.08\pm 0.28\) & \(81.12\pm 0.67\) & 6.2 \\ HAN (full batch)\({}^{*}\) & \(83.58\pm 1.46\) & \(94.04\pm 0.41\) & OOM & \(71.05\pm 2.26\) & OOM & \(83.27\pm 1.62\) & OOM & OOM & 9.6 \\ HAN (min batch) & \(80.77\pm 2.36\) & \(59.12\pm 0.52\) & \(26.05\pm 1.37\) & \(62.00\pm 0.46\) & \(85.71\pm 0.80\) & \(82.40\pm 2.68\) & \(63.0\pm 0.96\) & \(79.72\pm 0.62\) & 10.4 \\ AllDeepgsets & \(88.09\pm 1.52\) & \(59.68\pm 0.26\) & \(30.36\pm 1.67\) & \(67.82\pm 2.40\) & \(64.53\pm 0.33\) & \(80.70\pm 1.78\) & \(**78.46\pm 0.26**\) & \(81.06\pm 0.54\) & 5.6 \\ AllSetTransformer & \(88.69\pm 1.24\) & \(98.20\pm 0.20\) & \(\mathbf{36.89\pm 0.51}\) & \(69.33\pm 2.20\) & \(\mathbf{65.46\pm 0.25}\) & \(83.14\pm 1.92\) & \(\mathbf{78.46\pm 0.40}\) & \(81.38\pm 0.58\) & 3.1 \\ \hline PhenomNN\({}_{\text{simple}}\) & \(\mathbf{91.03\pm 1.04}\) & \(\mathbf{98.66\pm 0.20}\) & \(32.62\pm 0.40\) & \(\mathbf{71.77\pm 1.68}\) & \(\mathbf{64.41\pm 0.49}\) & \(\mathbf{86.96\pm 1.33}\) & \(\mathbf{78.46\pm 0.32}\) & \(\mathbf{81.74\pm 0.52}\) & **1.6** \\ PhenomNN & \(90.62\pm 1.88\) & \(98.61\pm 0.17\) & \(31.92\pm 0.36\) & \(70.71\pm 2.35\) & \(62.98\pm 1.36\) & \(85.28\pm 2.30\) & \(78.26\pm 0.26\) & \(81.41\pm 0.49\) & 3.1 \\ \hline \end{tabular}
\end{table}
Table 2: Results using the benchmarks from (Chien et al., 2022): Mean accuracy (%) \(\pm\) standard deviation. The number behind Walmart and House is the feature noise standard deviation for each dataset, and for HAN\({}^{*}\), additional preprocessing of each dataset is required (see (Chien et al., 2022) for more details). Boldfaced letters are used to indicate the best mean accuracy and underline is for the second. OOM indicates out-of-memory.
and the Major Key Project of PCL (PCL2021A12).
