Bayesian Sheaf Neural Networks

• Bayesian Sheaf Neural Networks are variational extensions of sheaf-based graph neural networks that treat the sheaf Laplacian as a latent random variable.
• They improve robustness on heterophilic graphs by replacing deterministic message passing with Bayesian averaging over learned sheaf restriction maps.
• Key techniques include variational inference using the ELBO, specialized distributions such as the Cayley transform for SO(d) maps, and robust optimization strategies.

Bayesian Sheaf Neural Networks (BSNNs) are variational Bayesian extensions of sheaf-based graph neural networks in which the sheaf Laplacian, or equivalently the restriction maps defining a cellular sheaf on a graph, is treated as latent random structure rather than as a deterministic point estimate. In this formulation, the network retains sheaf-diffusion message passing while replacing deterministic sheaf learning with variational inference over sheaves, with the stated aim of improving robustness on heterophilic graph data, especially when labeled data are limited (Gillespie et al., 2024). Within the broader sheaf-neural literature, deterministic Sheaf Neural Networks provide the mathematical substrate—stalks, restriction maps, global sections, and sheaf Laplacians—on top of which a Bayesian treatment can be built, but not every sheaf neural network is Bayesian in the technical sense (Mehrab et al., 29 Jan 2026).

1. Conceptual position within graph representation learning

BSNNs are motivated by the limitations of ordinary message passing on heterophilic graphs. Standard GNNs and many higher-order extensions work best on homophilic graphs, where neighboring nodes tend to share labels. On heterophilic graphs, adjacent nodes often belong to different classes, and ordinary message passing can oversmooth features and destroy class-separating information. Sheaf neural networks address this by replacing the graph Laplacian with a sheaf Laplacian $L_F$ or its normalized variant $\Delta_F$, so that neighboring nodes need not be aggregated through a single uniform averaging geometry (Gillespie et al., 2024).

The immediate precursor is deterministic Neural Sheaf Diffusion. In that setting, the sheaf is learned as a function of node features, which is expressive but also potentially overly sensitive to limited labeled data, random initialization, hyperparameter choices, and noise in the input features. BSNNs replace that point-estimate view by treating the sheaf itself as a latent random variable and learning a distribution over sheaves rather than a single sheaf. The resulting model is framed as a Bayesian counter-measure to brittle sheaf learning, not as a replacement for sheaf diffusion as such.

A recurring misconception is to equate any sheaf-based graph network with a Bayesian model. That is incorrect. Deterministic sheaf architectures may learn restriction maps, normalized sheaf Laplacians, and sheaf-aware diffusion operators without introducing priors, posteriors, variational objectives, or explicit uncertainty quantification. In the technical sense used by the BSNN literature, the defining Bayesian step is the introduction of a variational posterior $q_\phi(F\mid X,y)$ over sheaf structure.

2. Sheaf-theoretic and diffusion foundations

The sheaf-neural framework begins with an undirected graph $G=(V,E)$. In the deterministic SNN presentation, the graph is viewed as a finite topological space via the Alexandrov topology. Each vertex $v$ has an associated open star $U_v$, each edge $e$ has $U_e=\{e\}$, and a sheaf $F$ of vector spaces is specified by assigning a vector space $F(U_v)$ to each vertex, a vector space $\Delta_F$0 to each edge, and linear restriction maps

$\Delta_F$1

whenever $\Delta_F$2 is a vertex of $\Delta_F$3 (Mehrab et al., 29 Jan 2026).

If a graph with features is denoted $\Delta_F$4, each node $\Delta_F$5 carries a feature vector $\Delta_F$6. The collection $\Delta_F$7 can be regarded as sections of a constant sheaf $\Delta_F$8. The essential sheaf-theoretic departure from ordinary GNNs is that node-to-edge relations need not be identity maps: the restriction maps encode structured local compatibility. The usual sheaf principle then applies: local data can be “glued” into global data when the local pieces are consistent.

This consistency is measured through the sheaf Laplacian. With chain spaces

$\Delta_F$9

and a coboundary operator $q_\phi(F\mid X,y)$0 defined for an oriented edge $q_\phi(F\mid X,y)$1 connecting $q_\phi(F\mid X,y)$2 and $q_\phi(F\mid X,y)$3 by

$q_\phi(F\mid X,y)$4

the sheaf Laplacian is

$q_\phi(F\mid X,y)$5

It is positive semidefinite, and its kernel corresponds to global sections. When all stalks are $q_\phi(F\mid X,y)$6 and all restriction maps are identities, the construction reduces to the ordinary graph Laplacian $q_\phi(F\mid X,y)$7, so the sheaf model is a strict generalization of Laplacian-based GNNs (Mehrab et al., 29 Jan 2026).

The same structure underlies sheaf diffusion. In Neural Sheaf Diffusion, node features evolve under

$q_\phi(F\mid X,y)$8

where $q_\phi(F\mid X,y)$9. The harmonic space is

$G=(V,E)$0

and the diffusion limit satisfies

$G=(V,E)$1

Accordingly, sheaf diffusion learns a task-adapted local consistency geometry whose limiting space need not be the constant-signal subspace of the ordinary graph Laplacian (Dönmez et al., 11 May 2026).

3. Variational Bayesian formulation

The BSNN formulation is given for semi-supervised node classification on a graph $G=(V,E)$2 with node features $G=(V,E)$3, observed labels $G=(V,E)$4, and unobserved labels $G=(V,E)$5. The predictive model is $G=(V,E)$6, implemented by a sheaf neural network whose layers use sampled sheaves $G=(V,E)$7. The paper’s graphical model is summarized by

$G=(V,E)$8

where $G=(V,E)$9 is the adjacency matrix, $v$0 parameterizes the variational sheaf learner, and $v$1 parameterizes the classifier and sheaf-network weights (Gillespie et al., 2024).

The sheaf diffusion layer is based on the discretized process

$v$2

and the layer update used in the paper is

$v$3

with $v$4, $v$5, stalk dimension $v$6, feature-channel count $v$7, and nonlinearity $v$8. A stack of such layers is used, and each layer can have its own sampled sheaf $v$9.

The Bayesian step is to promote the sheaf to a latent variable. The intractable posterior $U_v$0 is approximated by a variational family $U_v$1. The evidence lower bound is written as

$U_v$2

with

$U_v$3

This is the objective optimized during training.

The variational sheaf learner is defined per incident node-edge pair $U_v$4. Both posterior and prior are factorized: $U_v$5 Their distributional parameters are produced by an MLP: $U_v$6 where $U_v$7 is the other endpoint of edge $U_v$8. This is an amortized posterior over restriction maps: the network infers local sheaf geometry from endpoint features rather than storing a single deterministic map for each incidence.

4. Restriction-map families and the Cayley-distribution construction

The paper treats three families of restriction maps. For general linear restriction maps, $U_v$9 is a $e$0 matrix, with $e$1. A latent Gaussian sample

$e$2

is reshaped into a $e$3 matrix, with reparameterization

$e$4

and prior

$e$5

For invertible diagonal restriction maps, $e$6, a Gaussian sample is drawn in $e$7, and the matrix is set to

$e$8

with standard normal prior on the diagonal vector.

The special orthogonal case is the paper’s distinctive probabilistic contribution. Here $e$9, $U_e=\{e\}$0 is a mean rotation, $U_e=\{e\}$1 is a concentration parameter, and the posterior family is

$U_e=\{e\}$2

with prior equal to the uniform distribution $U_e=\{e\}$3. The family $U_e=\{e\}$4 is defined through the Cayley transform

$U_e=\{e\}$5

whose inverse, for $U_e=\{e\}$6 with $U_e=\{e\}$7 not an eigenvalue, is

$U_e=\{e\}$8

This yields a reparameterizable probability distribution on $U_e=\{e\}$9, which is essential for gradient-based variational optimization (Gillespie et al., 2024).

The density with respect to normalized Haar measure is given in closed form: $F$0 For $F$1, the paper identifies $F$2 with a wrapped Cauchy distribution on the circle $F$3. For $F$4, it relates $F$5 to an angular central Gaussian distribution on $F$6, pushed forward through the standard two-fold cover $F$7.

The same paper also strengthens an expressivity claim for special orthogonal sheaves. It proves that for any $F$8, if $F$9 denotes the class of connected graphs with $F(U_v)$0 classes, then $F(U_v)$1 has linear separation power over $F(U_v)$2. This supports the use of $F(U_v)$3-valued restriction maps not only as a probabilistic convenience but also as an expressive sheaf family.

5. Optimization, inference, and empirical behavior

Training optimizes a Monte Carlo estimate of the ELBO. When the KL term is tractable, the objective is

$F(U_v)$4

and when it is not, the paper uses

$F(U_v)$5

A KL weight $F(U_v)$6 is cyclically annealed to prevent KL vanishing. For Gaussian posteriors against standard normal priors, analytic KL expressions are used; for Cayley versus uniform on $F(U_v)$7, analytic KL formulas are given for $F(U_v)$8, while for $F(U_v)$9 the KL is estimated by sampling (Gillespie et al., 2024).

The implementation pipeline is explicit. The variational sheaf learner first applies a linear layer to resize node features to dimension $\Delta_F$00, concatenates endpoints of each edge to build $\Delta_F$01 incident-pair inputs, uses an MLP to produce $\Delta_F$02, samples sheaf restriction maps $\Delta_F$03, propagates features through $\Delta_F$04 sheaf layers, applies a final linear classifier, and optimizes the ELBO estimator with backpropagation. At test time, prediction is performed by ensembling over multiple sampled sheaves so that uncertainty in the learned variational posterior is integrated out in the predictive point estimate.

The reported experiments are on the WebKB datasets Texas, Wisconsin, and Cornell. Node features are 1703-dimensional bag-of-words vectors, classes are five—student, project, course, staff, and faculty—and the graph sizes are 183 nodes and 325 edges for Texas, 183 nodes and 298 edges for Cornell, and 251 nodes and 515 edges for Wisconsin. The evaluation protocol uses 10 fixed splits of the data with $\Delta_F$05 train/validation/test proportions, hyperparameter grid search, retraining and evaluation with 30 random seeds for the best setting, and a Wilcoxon signed-rank test at $\Delta_F$06. The grid includes hidden channels $\Delta_F$07, stalk dimensions $\Delta_F$08, layers $\Delta_F$09, dropout $\Delta_F$10, learning rate $\Delta_F$11, ELU activation, weight decay $\Delta_F$12, patience $\Delta_F$13, max epochs $\Delta_F$14, Adam optimization, and ensemble size $\Delta_F$15 for BSNN.

Empirically, the paper reports that on each dataset the best test accuracy among the compared models is achieved by one of the BSNN variants; $\Delta_F$16-BSNN outperforms the corresponding deterministic $\Delta_F$17-NSD on all three datasets; Gen-BSNN significantly outperforms Gen-NSD on all three datasets; and diagonal results are mixed but still competitive. The standard deviation across 30 seeds is generally lower for BSNN than NSD, and when performance is averaged over all hyperparameter settings, BSNN consistently dominates the corresponding deterministic model. The stated interpretation is that Bayesian averaging over plausible sheaves is especially helpful when the sheaf family is expressive and training data are limited.

6. Deterministic sheaf neural networks and the non-Bayesian biomedical case

A separate sheaf-neural development presents Sheaf Neural Networks as a deterministic generalization of graph neural networks and applies them to a biomedical classification problem. That model uses a sheaf-aware diffusion operator in place of a standard graph convolution, with learned sheaves at each layer and general, unconstrained restriction maps. When the sheaf in every layer is the constant sheaf $\Delta_F$18, the architecture reduces to the ordinary Laplacian-based GNN case. Training is standard supervised node classification in a transductive setting, using Adam with a learning-rate scheduler for 400 epochs, no mini-batching, and regularization through dropout, batch normalization, layer normalization, and weight decay. Hyperparameters are selected by exhaustive grid search within each cross-validation fold, with the best setting chosen by mean validation accuracy across repeated 10-fold stratified cross-validation runs (Mehrab et al., 29 Jan 2026).

The biomedical application is osteosarcoma classification. The dataset contains 224 XANES spectra from bone tissue samples, each spectrum having 501 data points. The binary labels are 147 tumor samples and 77 control samples. Samples are reduced to 50 dimensions by PCA before graph construction; the graph is built by computing cosine similarity between every pair of samples in PCA space, sorting similarities, and adding edges from the most similar pairs until the graph becomes connected. Each node is therefore a patient tissue sample with a 50-dimensional feature vector, and the task is tumor-versus-healthy node classification.

Evaluation uses 10-fold stratified cross-validation repeated 5 times with different shuffles. The main metrics are fold accuracy and majority-vote accuracy; the appendix additionally reports precision, recall, F1, and AUC. The best-performing model is SheafGeneral, with fold accuracy

$\Delta_F$19

and majority-vote accuracy

$\Delta_F$20

For comparison, GraphSAGE achieves $\Delta_F$21 fold accuracy and $\Delta_F$22 majority-vote accuracy, while GAT and GCN are much lower, around $\Delta_F$23 and $\Delta_F$24, respectively. In the appendix, SheafGeneral also achieves precision $\Delta_F$25, recall $\Delta_F$26, F1 $\Delta_F$27, and AUC $\Delta_F$28.

This result is directly relevant to BSNNs because it demonstrates that the underlying deterministic sheaf architecture can outperform standard GNN baselines on a concrete biomedical task. It is not, however, a Bayesian sheaf neural network in the explicit sense of Bayesian inference, variational learning, priors, posteriors, or uncertainty quantification. The paper itself states that a Bayesian treatment could be built later on top of the deterministic sheaf-Laplacian framework.

7. Oversmoothing, representation geometry, and limitations

A later theoretical development reinterprets Neural Sheaf Diffusion through incidence-quiver representations. For an undirected graph $\Delta_F$29, the incidence quiver $\Delta_F$30 has vertices $\Delta_F$31 and arrows $\Delta_F$32 for every incidence $\Delta_F$33. Under the paper’s correspondence, finite-dimensional cellular sheaves on $\Delta_F$34 are equivalent to finite-dimensional representations of $\Delta_F$35, with fixed stalk dimensions giving a representation space

$\Delta_F$36

and gauge group

$\Delta_F$37

acting by base change

$\Delta_F$38

On this view, oversmoothing is not only spectral collapse but representation degeneration: direct-sum decompositions of the learned sheaf induce decompositions of the harmonic space

$\Delta_F$39

so collapse toward low-complexity summands can destroy discriminative information in the diffusion limit (Dönmez et al., 11 May 2026).

The same work studies the trivial subrepresentation $\Delta_F$40, in which all stalks align so that restriction maps act as identities on a one-dimensional subspace. On this summand, the Dirichlet energy reduces to ordinary graph Laplacian energy, and for a connected graph the corresponding global sections are constant. To bias training away from such collapse, the paper introduces the moment-map-inspired regularizers CentMM and ThetaMM: $\Delta_F$41

$\Delta_F$42

It also identifies an equal-stalk obstruction: when $\Delta_F$43, admissibility for learnable stability parameters forces the trivial all-object summand onto a stability wall, whereas non-uniform stalk dimensions such as $\Delta_F$44 remove this obstruction. On heterophilic benchmarks, the reported results are consistent with this mechanism: rectangular architectures can reduce variance or improve validation behavior, and on Wisconsin the rectangular ThetaMM model reaches $\Delta_F$45, while both rectangular and square models remain resilient up to $\Delta_F$46 before numerical instability appears at $\Delta_F$47 (Dönmez et al., 11 May 2026).

For Bayesian sheaf models, these results do not constitute a formal Bayesian derivation, but they sharpen the geometry of what a posterior over sheaves is implicitly ranging over. This suggests that priors over restriction maps should respect representation geometry rather than treat matrix entries as naive isotropic parameters; a plausible implication is that CentMM and ThetaMM can be read as soft geometric priors or MAP regularizers; and the equal-stalk obstruction suggests that the choice $\Delta_F$48 functions as a structural prior over feasible stability behavior, not merely as a capacity choice.

The current BSNN formulation also has explicit limitations. The variational posterior is factorized over incident node-edge pairs and therefore ignores posterior correlations between restriction maps; the Cayley family is exploited analytically most conveniently for $\Delta_F$49, while for $\Delta_F$50 KL terms are sampled rather than closed form; experiments are limited to WebKB node classification in a low-data regime; prediction cost is higher because inference uses ensembling over sampled sheaves; and the empirical study focuses on accuracy, standard deviation, and hyperparameter robustness rather than on calibration curves, predictive entropy, or out-of-distribution uncertainty analysis (Gillespie et al., 2024).